Ecuación del calor(Grupo Eau De Parfum(EDP)))

Trabajo realizado por estudiantes
Título	Ecuación del Calor. Grupo Eau De Parfum (EDP)
Asignatura	EDP
Curso	2023-24
Autores	Lestau Torres, Pablo López Rojo, Celia Muñoz Guijarro, Sofía
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

En este trabajo, exploramos el comportamiento térmico de una varilla metálica aislada lateralmente y sometida a condiciones de temperatura prescritas en sus extremos. Comenzamos planteando el sistema de ecuaciones diferenciales parciales que describe este fenómeno, considerando constantes de conducción térmica y calor específico unitarias para simplificar el análisis.

Luego, investigamos la solución del estado estacionario y su interpretación, así como la obtención de la solución mediante el método de separación de variables. Posteriormente, representamos gráficamente la evolución temporal de la temperatura en la varilla, observando cómo varía en diferentes instantes de tiempo y analizando el flujo de calor entrante y saliente en ambos extremos.

Exploramos también variaciones en el coeficiente de conductividad térmica y en las condiciones de contorno, investigando cómo estas modificaciones afectan la solución del problema. Finalmente, analizamos cómo cambia el comportamiento térmico de la varilla cuando se aísla uno de sus extremos y verificaremos el principio del máximo en el problema modificado.

Por último hemos explorado la utilidad de la solución fundamental en la resolución del problema de Cauchy asociado a la ecuación del calor.

2 Conocimientos previos

Antes de comenzar será de gran ayuda introducir una serie de conceptos que utilizaremos a lo largo este trabajo:

• Ley del calor de Newton: describe la tasa de pérdida o ganancia de calor de un objeto en relación con el ambiente circundante.La ley se expresa comúnmente como:

[math] \frac{d T(t)}{d t} = - k (T(t) - T_{\mathrm{a}}) [/math],

siendo [math]T(t) [/math]la temperatura del objeto en un tiempo dado, [math]T_{\mathrm{a}}[/math] la temperatura ambiente y [math]k[/math] la constante de proporcionalidad.

• Principio de conservación de la energía calorífica: establece que la variación de energía calorífica sobre un cuerpo V se debe al balance entre el calor que entra y sale del cuerpo más una producción externa.

• Ley de Fourier:intaura que el calor fluye desde regiones de alta temperatura a regiones de baja temperatura, y la cantidad de flujo de calor depende de la diferencia de temperatura y las propiedades del material.La expresión matemática de la Ley de Fourier es:

[math] \mathbf{q} = - k {\nabla} T [/math],

donde [math] \mathbf{q}[/math] es el vector de flujo de calor por unidad de superficie y [math] {\nabla} T [/math] es el gradiente del campo de temperatura en el interior del material.

• Coeficiente de conductividad térmica (D):cuantifica la capacidad de un material para conducir el calor. Se define como la cantidad de calor que pasa en un segundo a través de un metro cuadrado de material cuando la diferencia de temperatura a través del material,cuando la diferencia de temperatura a través del material es de un grado(ºC) es de un grado, y el espesor del material es de un metro.
• Principio del máximo: Sea [math]u \in C^{2,1}(Q_T)\cap (\overline{Q_T})[/math] tal que [math]u_t - \Delta u \leq 0[/math] en [math]Q_T[/math]. Entonces u alcanza su máximo en la frontera parabólica [math]\partial _p Q_T[/math]:

[math]\max_{\overline{Q_T}} u = \max_{\partial _p Q_T} u[/math]

3 Problema original

Habiendo ya introducido los conocimientos básicos, comenzaremos con nuestro problema. Consideramos una varilla metálica de 1 metro de longitud , que se encuentra aislada por su superficie lateral, de manera que la conducción de calor únicamente se produce en la dirección longitudinal. La temperatura inicial de la varilla es de 0ºC. En el extremo izquierdo se consigue mantener la temperatura a 0 ºC mientras que en el derecho la temperatura es siempre de 1 ºC. Además el calor específico toma un valor constante de 1.En cuanto al coeficiente de conductividad térmica,vamos a observar cómo sería nuestro problema para un valor cualquiera y seguidamente especificaremos el valor constante de 1.

Pero antes de abordar matemáticamente nuestro problema, vamos a intentar deducir lo que podría ocurrir. Considerando que la temperatura en el extremo derecho de la varilla es superior a la del extremo izquierdo, el calor comenzará a transferirse desde la zona más cálida, elevando la temperatura a lo largo de la varilla y promoviendo el desplazamiento del calor hacia la región más fría. Simultáneamente, el incremento de la temperatura interna provoca una reducción del flujo de calor entrante y a un aumento del saliente con el tiempo. Anticipamos que, en algún momento, ambos flujos se equilibrarán, conduciendo la temperatura de la varilla a un estado estacionario.

Por lo que nuestro problema a resolver queda:

[math] \left\{ \begin{aligned} &u_t - Du_{xx} = 0, & 0 \lt x \lt 1, t \gt 0, \\ &u(0, t) = 0, & t \gt 0, \\ &u(1, t) = 1, & t \gt 0, \\ &u(x, 0) = 0, & 0 \lt x \lt 1, \end{aligned} \right. [/math]

Debemos reescalar en tiempo, para facilitarnos los cálculos, por lo que realizamos el cambio de variable [math]\tau = \alpha t, \tau \in\mathbb{R} [/math]. Por tanto, [math] u(x,t)=U(x,\tau)[/math] y entonces [math] \alpha U_{\tau}-DU_{xx}=0[/math] .Tomando [math] \alpha = D[/math] , obtenemos el problema

[math] \left\{ \begin{aligned} &U_{\tau} - U_{xx} = 0, & 0 \lt x \lt 1, \tau \gt 0, \\ &U(0, \tau) = 0, & \tau \gt 0, \\ &U(1, \tau) = 1, & \tau \gt 0, \\ &U(x, 0) = 0, & 0 \lt x \lt 1, \end{aligned} \right. [/math]

Este problema se resuelve mediante el método de variación de las constantes, pero un factor clave para su implementación es que las condiciones de contorno estén homogeneizadas, es decir, que sean igual a cero. Esto se conseguirá mediante un cambio de variable en el que es imprescindible la llamada solución estacionaria. La solución estacionaria de nuestro problema es aquella que no varía con el tiempo [math] v=\lim_{\tau \to \infty} U(x,\tau)[/math] y por tanto, su derivada respecto al tiempo es nula. Nuestro sistema queda ahora:

[math] \left\{ \begin{aligned} &-v_{xx} = 0, & 0 \lt x \lt 1, \\ &v(0)= 0 \\ &v(1) = 1, \end{aligned} \right. [/math]

Resolviendo obtenemos que [math] v(x)=x [/math] es nuestra solución estacionaria, cuya expresión gráfica es:

import matplotlib.pyplot as plt
import numpy as np
# Define el rango de x entre 0 y 1
x = np.linspace(0, 1, 100)  
# Define la función  estacionaria v(x) = x
v_x = x
# Crea la gráfica
plt.plot(x, v_x)
# Etiquetas de los ejes
plt.xlabel('x')
plt.ylabel('v(x)')
# Título de la gráfica
plt.title('Gráfica de la solución estacionaria v(x)=x')
# Muestra la gráfica
plt.show()

Una vez hallada, podemos homogeneizar nuestro problema mediante el cambio de variable [math]w(x,t)= U(x,\tau)-v(x)[/math].Por tanto, sustituyendo en el sistema inicial :

[math] \left\{ \begin{aligned} &w_{\tau} - w_{xx} = 0, & 0 \lt x \lt 1, \tau \gt 0, \\ &w(0, \tau) = 0, & \tau \gt 0, \\ &w(1, \tau) = 1-1=0, & \tau \gt 0, \\ &w(x, 0) = 0-x=-x, & 0 \lt x \lt 1, \end{aligned} \right. [/math]

Nuestro sistema es ya homogéneo,por lo que aplicamos el método de separación de variables (que consiste en expresar la solución como el producto de una función que depende únicamente de la variable espacial x y una función dependiente de la variable temporal) y obtenemos como solución homogénea:

[math]w(x,\tau) = \sum_{k=1}^{\infty} \frac{2(-1)^{k}}{k\pi} e^ {-k^2\pi ^2 \tau}\sin(k\pi x)[/math]

Y deshaciendo los dos cambios de variable realizados, nuestra solución al problema original es por tanto:

[math]u(x,t) = x+ \sum_{k=1}^{\infty} \frac{2(-1)^{k}}{k\pi} e^ {-k^2\pi ^2 Dt }\sin(k\pi x)[/math]

Como en nuestro caso, el coeficiente de conductividad térmica es igual a 1, nuestro sistema es:

[math] \left\{ \begin{aligned} &u_t - u_{xx} = 0, & 0 \lt x \lt 1, t \gt 0, \\ &u(0, t) = 0, & t \gt 0, \\ &u(1, t) = 1, & t \gt 0, \\ &u(x, 0) = 0, & 0 \lt x \lt 1, \end{aligned} \right. [/math]

y nuestra solución:

[math]u(x,t) = x+ \sum_{k=1}^{\infty} \frac{2(-1)^{k}}{k\pi} e^ {-k^2\pi ^2 t }\sin(k\pi x)[/math]

Mostramos gráficamente la solución, observando que se cumplen las condiciones de frontera. Además según avanza el tiempo, la solución deja de depender de este, alcanzando un estado estacionario determinado justamente por la solución estacionaria [math] v(x)=x [/math]. Cabe destacar que en [math] t=0 [/math] obtenemos las oscilaciones correspondientes a la serie de Fourier de la condición inicial, que se aproxima mediante senos.

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

# Parámetros
L = 1       # Longitud de la varilla
T = 1      # Tiempo final
D=1       # Coeficiente de conductividad térmica
N = 10      # Número de términos en la serie
x_points = 100  # Puntos en la dirección x
t_points = 100  # Puntos en la dirección t

# Espacios
x = np.linspace(0, L, x_points)
t = np.linspace(0, T, t_points)
X, T = np.meshgrid(x, t)

# Solución usando los primeros N términos de la serie
def u(x, t, N):
    sum = 0
    for k in range(1, N+1):
        lambda_k = k**2 * np.pi**2 
        #Calculasmos los coeficientes C_k
        C_k = 2 * (-1)**(k) / (k * np.pi)  
        sum += C_k * np.sin(k*np.pi * x) * np.exp(-lambda_k * t*D) 
    return x+sum

# Calcular la solución en la malla
Z = u(X, T, N)

# Graficar
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111, projection='3d')
surf=ax.plot_surface(X, T, Z, cmap='viridis')

ax.set_xlabel('Posición x (m)')
ax.set_ylabel('Tiempo t (s)')
ax.set_zlabel('Temperatura u(x,t) (ºC)')
ax.set_title('Evolución de la Temperatura en la Varilla')

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

3.1 Flujo

Para estudiar el flujo de calor saliente y entrante en ambos extremos de la varilla a lo largo del tiempo, nos debemos basar en la ley de Fourier [math] \mathbf{q} = - k u_x(x,t) [/math] Por lo que necesitamos obtener la derivada de la solución con respecto al espacio:

[math] u_x(x,t)= 1+\sum_{k=1}^{\infty}(-1)^{k+1}2 e^{-k^2\pi ^2 t }\cos(k\pi x) [/math]

Ahora, el flujo en los extremos es:

[math] u_x(0,t)=1 +\sum_{k=1}^{\infty}(-1)^{k+1}2 e^{-k^2\pi ^2 t } [/math]

[math] u_x(1,t)=1-\sum_{k=1}^{\infty}2 e^{-k^2\pi ^2 t } [/math]

Gráficamente,

import numpy as np
import matplotlib.pyplot as plt
# Parámetros
L = 1       # Longitud de la varilla
T = 1      # Tiempo final
D=1       # Coeficiente de conductividad térmica
N = 10      # Número de términos en la serie
x_points = 100  # Puntos en la dirección x
t_points = 100  # Puntos en la dirección t

# Espacios
x = np.linspace(0, L, x_points)
t = np.linspace(0, T, t_points)
X, T = np.meshgrid(x, t)

# Derivada de nuestra solución respecto al espacio, multiplicado por -k
def u(x, t, N):
    sum = 0
    for k in range(1, N+1):
        lambda_k = k**2 * np.pi**2 
        #Calculamos los coeficientes C_k mediante extensión impar
        C_k = 2 * (-1)**(k+1)  
        sum += C_k * np.cos(k*np.pi * x) * np.exp(-lambda_k * t) 
    return 1+sum

# Calcular la solución en la malla
Z = u(X, T, N)
# Graficar
plt.figure(figsize=(8, 6))
plt.plot(t, u(0,t,N), label='Flujo en x=0', color='blue')  # Cálculo del flujo en x=0
plt.plot(t, u(1,t,N), label='Flujo en x=1', color='red')   # Cálculo del flujo en x=1
plt.xlabel('Tiempo [s]')
plt.ylabel(r'Flujo de Calor [$\frac{^\circ C}{m^2 \cdot s}$]')
plt.title('Flujo de Calor en los Extremos de la Varilla ')
plt.legend()
plt.grid(True)
plt.show()

Al comenzar con un flujo de calor nulo en el extremo izquierdo, esto indica que hay una gran diferencia de temperatura entre ambos extremos de la varilla. El extremo izquierdo está perdiendo calor rápidamente hacia la zona más fría. A medida que el calor se va transfiriendo por la varilla, el gradiente de la temperatura cambia y el extremo derecho comienza a calentarse, debido al flujo de calor que proviene del otro extremo. Este cambio conduce a que el flujo en ambos extremos comience a estar en equilibrio. Esta estabilización de flujo en -1 refleja que se ha alcanzado un estado estacionario en el sistema, indicando que cualquier calor entrante en la varilla, se compensa con el calor saliente, manteniendo una distribución de temperatura uniforme.

Este comportamiento refleja la Ley de Fourier, donde al principio la gran diferencia de temperatura induce un flujo de calor bastante significativo, que se va estabilizando llegando a un flujo de -1 en ambos extremos, alcanzando así un gradiente de temperatura constante, con la misma cantidad de calor saliendo hacia fuera desde ambos extremos.

4 Cambio de coeficiente de conductividad térmica

Hemos estudiado el problema de Cauchy para un coeficiente de difusión de 1,veamos que ocurre si lo disminuimos a su mitad, es decir,[math]D=\frac{1}{2}[/math].Nuestro nuevo problema por tanto es:

[math] \left\{ \begin{aligned} &u_t - \frac{1}{2}u_{xx} = 0, & 0 \lt x \lt 1, t \gt 0, \\ &u(0, t) = 0, & t \gt 0, \\ &u(1, t) = 1, & t \gt 0, \\ &u(x, 0) = 0, & 0 \lt x \lt 1, \end{aligned} \right. [/math]

Cuya solución es

[math]u(x,t) = x+ \sum_{k=1}^{\infty} \frac{2(-1)^{k}}{k\pi} e^ {\frac{-k^2\pi ^2 t}{2}}\sin(k\pi x)[/math]

Gráficamente, únicamente hemos cambiado el valor de D en el código (2) y obtenemos:

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# Parámetros
L = 1       # Longitud de la varilla
T = 1      # Tiempo final
D=1/2       # Coeficiente de conductividad térmica
N = 10      # Número de términos en la serie
x_points = 100  # Puntos en la dirección x
t_points = 100  # Puntos en la dirección t
# Espacios
x = np.linspace(0, L, x_points)
t = np.linspace(0, T, t_points)
X, T = np.meshgrid(x, t)
# Solución usando los primeros N términos de la serie
def u(x, t, N):
    sum = 0
    for k in range(1, N+1):
        lambda_k = k**2 * np.pi**2 
        #Calculasmos los coeficientes C_k
        C_k = 2 * (-1)**(k) / (k * np.pi)  
        sum += C_k * np.sin(k*np.pi * x) * np.exp(-lambda_k * t*D) 
    return x+sum
# Calcular la solución en la malla
Z = u(X, T, N)
# Graficar
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111, projection='3d')
surf=ax.plot_surface(X, T, Z, cmap='viridis')
ax.set_xlabel('Posición x (m)')
ax.set_ylabel('Tiempo t (s)')
ax.set_zlabel('Temperatura u(x,t) (ºC)')
ax.set_title('Evolución de la Temperatura en la Varilla')
fig.colorbar(surf, shrink=0.5, aspect=5)
plt.show()

No se puede apreciar ninguna diferencia aparente entre las soluciones para los diferentes coeficientes. Para observar algún cambio vamos a dibujar la diferencia entre la solución y el estado estacionario en [math] x = 1/2[/math] para [math] t \in [0, 1][/math], es decir [math]u(1/2, t) − v(1/2)[/math],

import numpy as np
import matplotlib.pyplot as plt
x=0.5
# Definir la función de temperatura con un número limitado de términos en la serie
def u(x, t, alpha, n_terms=10):   
    sum_series = sum(2 * (-1)**(k + 1) / (k * np.pi) * np.sin(k * np.pi * x) * np.exp(-k**2 * np.pi**2 * alpha * t) for k in range(1, n_terms + 1))
    return  x - sum_series
# Intervalo de tiempo
t = np.linspace(0, 1, 100)
# Calcular la temperatura en x=1/2 para diferentes alphas
u_k1 = u(0.5, t, alpha=1)-x  # D=1, correspondiente a alpha=1
u_k_half = u(0.5, t, alpha=0.5)-x  # D=1/2, correspondiente a alpha=0.5
# Graficar
plt.plot(t, u_k1, label='k=1', color='blue')
plt.plot(t, u_k_half, label='k=1/2', color='red')
plt.xlabel('Tiempo (s)')
plt.ylabel('u(1/2, t) − v(1/2)')
plt.title('Efecto del Coeficiente de Conductividad Térmica ')
plt.legend()
plt.show()

El coeficiente de conductividad térmica juega un papel fundamental en determinar la velocidad a la que se propaga la temperatura dentro de un material. Al disminuir este coeficiente a la mitad, el proceso de conducción térmica a través de la varilla experimenta una notable ralentización. Esto se traduce en que, en cualquier instante dado t, la distribución de la temperatura a lo largo de la varilla VA avanzando a la mitad de velocidad en comparación con un escenario donde el coeficiente de conductividad térmica es de 1. Dicho cambio incide directamente en la rapidez con que se logra el equilibrio térmico en la varilla, Resultando en una difusión térmica más lenta hacia el equilibrio.

5 Cambio de condiciones

Supongamos ahora una barra unidimensional con extremos fijos a una temperatura de [math]0^\circ C[/math]. La temperatura inicial en cualquier punto de la barra está dada por la función [math]u(x,0) = max\{0, 1−4|x−1/2|\}[/math], lo que significa que inicialmente la temperatura en la barra varía de acuerdo con esta función.

[math] \left\{ \begin{aligned} &u_t - u_{xx} = 0, & 0 \lt x \lt 1, t \gt 0, \\ &u(0, t) = 0, & t \gt 0, \\ &u(1, t) = 0, & t \gt 0, \\ &u(x, 0) = max\{0, 1−4|x−1/2|\}, & 0 \lt x \lt 1, \end{aligned} \right. [/math]

Supongamos que, en tiempo infinito, \( t \rightarrow \infty \), se llega a una solución estacionaria y que, por tanto, ya no varía en tiempo. De modo, que tomando esta suposición llegamos al siguiente sistema con \( u(x,t) \xrightarrow{\text{t tiende a }\infty} v(x) \).

[math] \left\{ \begin{aligned} &v_{xx} = 0, & 0 \lt x \lt 1, \\ &v(0) = 0, & \\ &v(1) = 0, & \end{aligned} \right. [/math]

Cuya solución es la constante \( v(x)=0 \). De modo que el sistema homogeneizado, coincide con el original, pues \( w(x,t)=u(x,t)-v(x) =u(x,t)-0=u(x,t)\) Antes de exigir la condición inicial, la solución coincide con la del primer apartado, pues nos encontramos ante la misma ecuación diferencial, salvo por la función \( u(x,0)\), la cual no toma parte hasta el final de la resolución de esta suponiendo separación de variables. De modo que

[math]u(x,t) = \sum_{k=1}^{\infty} C_k e^{(-k^2\pi ^2 t )}\sin(k \pi x)[/math]

Para hallar el valor de los coeficientes [math] C_k[/math], buscamos la serie de Fourier de la función [math] max\{0, 1−4|x−1/2|\}[/math] extendida de forma impar. Esto es posible, ya que

[math] \int^1_0 (max\{0, 1−4|x−1/2|\})^2 dx=0.1667 \lt \infty[/math]

El cálculo de los coeficientes se ha realizado numéricamente aproximando por la fórmula del trapecio. Los coeficientes se pueden consultar ejecutando el siguiente código en Python

import numpy as np 
def f(x):   # Defino función que se quiere aproximar
    return np.maximum(0, 1 - 4 * np.abs(x - 1/2))

def fourier_coefficients(f, n, x):
    coefficients = np.zeros(n) # Guardo los coeficientes en una lista
    for k in range(1,n+1):     # Compiezo en k=1, pues sin(0)=0
        coefficients[k-1] = 2 * np.trapz(f(x) * np.sin(k * np.pi * x), x)
    return coefficients

fourier_coefficients(f, 10, np.linspace(0, 1, 1000))

De modo que ya se ha calculado la solución de este nuevo problema. Visualmente, responde a la siguiente gráfica, resultado de ejecutar este código:

import sympy as sp
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Parámetros
L = 1       # Longitud de la varilla
T = 1       # Tiempo final
N = 10      # Número de términos en la serie
# Espacios
x = np.linspace(0, L, 1000) # Puntos en la dirección x
t = np.linspace(0, T, 1000) # Puntos en la dirección t
X, T = np.meshgrid(x, t)    # Generar una cuadrícula de puntos

def f(x):    # Defino mi función condición inicial
    return np.maximum(0, 1 - 4 * np.abs(x - 1/2))

def fourier_coefficients(f, n, x): # Coeficientes, recuperado del anterior trabajo
    coefficients = np.zeros(n)
    for k in range(1,n+1):         # Comenzamos en 1, pues sen(0)=0
        coefficients[k-1] = 2 * np.trapz(f(x) * np.sin(k * np.pi * x), x)
    return coefficients.tolist()   # Convertir a lista
coef=fourier_coefficients(f, N, x)

def u(x, t, n): # Definir la solución en términos de la serie de Fourier con N términos
    sum = 0                 # Almacenamos la suma (resultado) hasta el término 10
    for k in range(1, n+1): # Para cada k calculo el lambda y sumo su coeficiente respectivo
        lambda_k = k* np.pi
        sum += coef[k-1] * np.sin(lambda_k * x) * np.exp(-lambda_k**2 * t)
    return sum

Z = u(X, T, N) # Calcular los valores de la función para cada par (x, t)

# Graficar la función en 3D
fig = plt.figure(figsize=(12, 8))  # Cambio aquí para hacer la figura más grande
ax = fig.add_subplot(111, projection='3d')
surf=ax.plot_surface(X, T, Z,cmap='viridis')

ax.set_xlabel('Posición x [m]')
ax.set_ylabel('Tiempo t [s]')
ax.set_zlabel('Temperatura u(x,t) [ºC]')
ax.set_title('Evolución de la Temperatura en la Varilla')
fig.colorbar(surf, shrink=0.5, aspect=5)  #Añado una guía de colores
surf.set_clim(0, 1)                       # Le introduzco cota inf/sup a los colores
plt.show()

​ Vemos como la temperatura sufre una gran variación en tiempos cercanos a 0, y rápidamente tiende a su estado estacionario, [math]v(x)=0[/math]. Para poder visualizarlo mejor mostraremos una un gif donde la imagen cambia a medida que el tiempo pasa, restringido a [math]t \in [0,0.3][/math], pues a partir de ahí, la temperatura está suficientemente cerca de la estacionaria y, por tanto, no hay gran variación en la temperatura.

En un inicio, en \(t=0\) la temperatura toma valores negativos. Esto debido a que la exponencial se anula y son los coeficientes de Fourier quienes exigen el signo de la función. Una vez que [math]t\gt0[/math], u se hace positivo en todas partes, lo que significa que la temperatura se difunde casi instantáneamente a lo largo del eje x, resultando en una velocidad de propagación infinita. Es decir, cualquier perturbación se siente enseguida en toda parte del dominio.

Terminamos el estudio de esta ecuación de difusión con el estudio del flujo, la inversa de la variación de la variable x de la solución. Para ello derivamos la solución que hemos obtenido. En términos de [math]C_k[/math], puesto que no son dependientes de x, obtenemos: [math]\frac{\partial}{\partial x}u(x,t)=\sum_{k=1}^{\infty} C_k e^{(-k^2\pi ^2 t )}\cos(k \pi x)\frac{1}{k \pi}[/math] como flujo de la ecuación homogénea. Adjuntamos el código con el cual se ha calculado la expresión de forma numérica y su respectiva representación gráfica:

import sympy as sp
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Parámetros
L = 1       # Longitud de la varilla
T = 1       # Tiempo final
N = 10      # Número de términos en la serie
# Espacios
x = np.linspace(0, L, 1000) # Puntos en la dirección x
t = np.linspace(0, T, 1000) # Puntos en la dirección t
X, T = np.meshgrid(x, t)

def f(x):   # Defino la función de la ci para hallar los coef
    return np.maximum(0, 1 - 4 * np.abs(x - 1/2))

def fourier_coefficients(f, n, x): # Hallo los coef de Fourier  
    coefficients = np.zeros(n)
    for k in range(1,n+1): 
        coefficients[k-1] = 2 * np.trapz(f(x) * np.sin(k * np.pi * x), x)
    return coefficients.tolist()   # Convertir a lista
coef=fourier_coefficients(f, N, x)

def u(x, t, n):  # Definir la derivada en términos de la serie de Fourier con N términos
    sum = 0
    for k in range(1, n+1):
        lambda_k = k* np.pi
        sum += coef[k-1] *np.cos(k * np.pi * x)/(k * np.pi) * np.exp(-lambda_k**2 * t)
    return sum

# Graficamos
plt.figure(figsize=(8, 6))
plt.plot(t, u(0,t,N), label='Flujo en x=0', color='blue')  # Cálculo del flujo en x=0
plt.plot(t, u(1,t,N), label='Flujo en x=1', color='red')   # Cálculo del flujo en x=1
plt.xlabel('Tiempo [s]')
plt.ylabel(r'Flujo de Calor [$\frac{^\circ C}{m^2 \cdot s}$]')
plt.title('Flujo de Calor en los Extremos de la Varilla ')
plt.legend()
plt.grid(True)
plt.show()

Observamos dos casos parecidos: Uno desciende hasta alcanzar el 0 (rojo) y otro asciende (azul). El primer caso, un flujo positivo que disminuye a cero, significa que el calor está siendo transferido a lo largo de la varilla en una dirección específica, pero la velocidad de transferencia de calor disminuye gradualmente hasta que alcanza cero, lo que indica un equilibrio térmico en la varilla. En el segundo caso, el flujo de calor comienza siendo negativo y alcanza cero, esto implicaría que al principio está ocurriendo una pérdida de calor en la dirección específica a lo largo de la varilla, pero esta pérdida disminuye hasta que se detiene completamente.

6 Varilla con extremo aislado

En este estudio modelamos la difusión del calor en un alambre. Inicialmente, consideramos un alambre con una temperatura inicial de \( 0^\circ C \) en su extremo izquierdo. En contraste, el extremo derecho del alambre está aislado térmicamente, lo que significa que no hay transferencia de calor a través de este. En nuestro análisis, introducimos una fuente de calor en el centro del alambre en el instante inicial. Luego, examinamos cómo se propaga y se distribuye esta fuente de calor a lo largo del alambre con el tiempo.

[math] \left\{ \begin{aligned} &u_t - u_{xx} = 0, & 0 \lt x \lt 1, t \gt 0, \\ &u(0, t) = 0, & t \gt 0, \\ &u_x(1, t) = 0, & t \gt 0, \\ &u(x, 0) = max\{0, 1−4|x−1/2|\}, & 0 \lt x \lt 1, \end{aligned} \right. [/math]

Que la barra se encuentra aislada térmicamente en su extremo en \( x = 1 \) implica que no hay flujo de calor a través de este extremo. La ley de Fick, establece que el flujo de calor es proporcional al gradiente de temperatura según la ecuación \( q = -k \frac{du}{dx} \), donde \( q \) es el flujo de calor, \( k \) es la conductividad térmica y \( \frac{du}{dx} \) es el gradiente de temperatura a lo largo de la barra. En este contexto, si el flujo en dicho extremo es nulo, el gradiente de temperatura es nulo, por tanto \( u_x(1, t) = 0 \). Observamos además que la solución estacionaria en nuestro caso es \( v(x) = 0 \), pues es la solución del problema de contorno.

[math] v''(x) = 0, \quad v(0) = 0, \quad v'(1) = 0 [/math]

De forma análoga a como ocurría en el apartado anterior, el problema que estamos estudiando en este caso es homogéneo. Por tanto, sin necesidad de hacer ningún cambio, procedemos a resolverlo mediante el método de separación de variables. Suponiendo que \( u(x, t) = X(x) \cdot T(t) \). Obtenemos como solución \( X(x) = A \sin(\pi(1/2+k)x) \) y \( T(t) = B e^{-(\pi^2)(1/2+k)^2 t} \). De modo que, la solución sin tener en cuenta la condición inicial es:

[math] u_k(x, t) = c_k \cdot \sin(\pi(1/2+k)x) \cdot e^{-(\pi^2)(1/2+k)^2 t} [/math]

Para asegurar la existencia de solución del problema, es necesario verificar si existen ciertos coeficientes \( c_k \) tal que la solución pueda ser definida como:

[math] u(x, t) = \sum_{k=0}^{\infty} c_k \cdot \sin(\pi(1/2+k)x) \cdot e^{-(\pi^2)(1/2+k)^2 t} [/math]

Para ello, es necesario verificar que las autofunciones \( \{\sin(\pi(1/2+k)x)\}_{k=0}^{\infty} \) forman un conjunto completo en \( L^2 [0, 1]\) y ortogonal.

Para comprobar la ortogonalidad entre las autofunciones, calculamos la integral de su producto, o lo que es igual a calcular el producto escalar de \( L^2 [0,1]\), para dos funciones distintas \( \sin(\pi(1/2+k)x) \) y \( \sin(\pi(1/2+m)x) \) con \( k \neq m \in \mathbb{N}\):

[math] \int_{0}^{1} \sin(\pi(1/2+k)x) \sin(\pi(1/2+m)x) \, dx = \frac{\sin((k - m) \pi)}{2 (k - m) \pi} - \frac{\sin((1 + k + m) \pi)}{2 (1 + k + m) \pi} [/math]

Esta integral se anula para \( k,m \in \mathbb{N} \), lo que demuestra que las autofunciones son ortogonales entre sí en el intervalo \([0, 1]\). Por otro lado, dado que demostrar que las funciones forman una base completa para el espacio \( L^2[0,1] \) puede ser un proceso complejo y requiere un análisis detallado, en este trabajo supondremos que la forman.

import numpy as np
import matplotlib.pyplot as plt

# Parámetros
L = 1       # Longitud de la varilla
T = 1       # Tiempo final
N = 10      # Número de términos en la serie
# Espacios
x = np.linspace(0, L, 1000) # Puntos en la dirección x
t = np.linspace(0, T, 1000) # Puntos en la dirección t
X, T = np.meshgrid(x, t)    # Generar una cuadrícula de puntos

def f(x):   # Defino mi fucnión condición inicial
    return np.maximum(0, 1 - 4 * np.abs(x - 1/2))


def fourier_coefficient(f, k, x):  # Coeficientes, recuperado del anterior trabajo
    ck = 2* np.trapz(f(x) * np.sin(np.pi*(1/2+k)*x) ,x)
    return ck

# Definir la solución en términos de la serie de Fourier con N términos
def u(x, t, N):
    sum = 0                 # Almacenamos la suma hasta el término 10
    for k in range(0, N+1): # Para cada k cálculo su coeficiente y lo sumo al resultado
        ck = fourier_coefficient(f, k, x)
        sum += ck * np.sin(np.pi*(1/2+k)*x) * np.exp(-(np.pi**2)*(1/2+k)**2 * t)
    return sum


Z = u(X, T, N)  # Calcular los valores de la función para cada par (x, t)

# Graficar la función en 3D
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111, projection='3d')
surf=ax.plot_surface(X, T, Z,cmap='viridis')

ax.set_xlabel('Posición x [m]')
ax.set_ylabel('Tiempo t [s]')
ax.set_zlabel('Temperatura u(x,t) [ºC]')
ax.set_title('Evolución de la Temperatura en la Varilla')
fig.colorbar(surf, shrink=0.5, aspect=5)  #Añado una guía de colores
surf.set_clim(0, 1)                       # Le introduzco cota inf/sup a los colores
plt.show()

Gracias a los desarrollos en serie de Fourier, podemos calcular los coeficientes \( c_k \) para cumplir así la condición inicial. Dichos coeficientes se calculan mediante la integral:

[math] c_k := \frac{\int_{0}^{1} \max\{0, 1 - 4(x - 1/2)\} \sin(\pi(1/2+k)x) \, dx}{\int_{0}^{1} \sin^2(\pi(1/2+k)x) \, dx} [/math]

Dado que la integral es difícil de calcular analíticamente y en consecuencia es difícil obtener una expresión dependiente del parámetro [math] k[/math], la resolveremos de forma numérica mediante el método del trapecio, de forma análoga a como se hizo en el modelo sin aislante térmico previo. Para una comprensión más completa de la evolución de la temperatura del alambre en cuestión, es importante realizar un análisis visual de la solución propuesta. La representación gráfica de la solución en el intervalo de tiempo \( t \in [0, 1] \), considerando los primeros 10 términos de la serie de Fourier, ofrece una visualización tridimensional de la distribución del calor a lo largo del alambre en función del tiempo.

Es relevante destacar que a medida que transcurre el tiempo, la solución converge gradualmente hacia su estado estacionario, es decir, \( v(x)=0 \) . Además, observamos que se verifica el principio del máximo, dado que tanto la temperatura máxima y mínima se alcanza en la frontera parabólica. Esta última afirmación la podemos verificar de forma sencilla, obteniendo ambos valores de la función. El código utilizado mostró los siguientes datos.

# Calcular la temperatura máxima y mínima
max_temperature = np.max(Z)
min_temperature = np.min(Z)


print("Temperatura máxima:", max_temperature)
print("Temperatura mínima:", min_temperature)


# Encontrar las coordenadas donde se alcanza la temperatura máxima y mínima
max_index = np.unravel_index(np.argmax(Z), Z.shape)
min_index = np.unravel_index(np.argmin(Z), Z.shape)


max_x, max_t = x_values[max_index[1]], t_values[max_index[0]]
min_x, min_t = x_values[min_index[1]], t_values[min_index[0]]


print("Coordenadas de la temperatura máxima: (x={}, t={})".format(max_x, max_t))
print("Coordenadas de la temperatura mínima: (x={}, t={})".format(min_x, min_t))

Resultados de temperatura máxima y mínima

Parámetro	Valor
Temperatura máxima	\(0.92 ^\circ \text{C}\)
Temperatura mínima	\(-0.03 ^\circ \text{C}\)

Coordenadas de temperatura máxima y mínima

Coordenadas	Valor
Temperatura máxima	\((x=0.50, t=0.00)\)
Temperatura mínima	\((x=0.81, t=0.00)\)

Observamos que el máximo debería ser de [math] 1^\circ \text{C} [/math]. No obstante, el resultado reportado por el código es [math] 0.92 ^\circ \text{C} [/math]. Esto se debe al uso de la aproximación por series de Fourier. Lo mismo ocurre con el valor mínimo que debería se de [math] 0^\circ \text{C} [/math], no obstante el código reporta [math]-0.03^\circ \text{C} [/math]. Además dichos puntos están en la frontera parabólica, como afirma el principio del máximo.

En este punto es interesante verificar si nuestro modelo es correcto. En tal caso, el extremo izquierdo de la barra no debería perder calor. Por lo tanto, calculamos la derivada de nuestra expresión y, gracias a la Ley de Fick, sabemos que el flujo de calor es igual a la derivada cambiada de signo. Obtenemos los siguientes resultados:

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

import numpy as np
import matplotlib.pyplot as plt

# Parámetros
L = 1       # Longitud de la varilla
T = 1       # Tiempo final
N = 100     # Número de términos en la serie
# Espacios
x_values = np.linspace(0, L, 1000) # Puntos en la dirección x
t_values = np.linspace(0, T, 1000) # Puntos en la dirección t

def f(x):   # Defino mi función condición inicial
    return np.maximum(0, 1 - 4 * np.abs(x - 1/2))

def fourier_coefficient(f, k, x):   # Coeficientes, recuperado del anterior trabajo
    ck = 2* np.trapz(f(x) * np.sin(np.pi*(1/2+k)*x) ,x)
    return ck

def derivative_solution_x(x, t, N): # Definir la derivada de la solución respecto a x
    result = 0
    for k in range(0, N+1):
        ck = fourier_coefficient(f, k, x_values)
        result += ck * np.pi * (1/2 + k) * np.cos(np.pi * (1/2 + k) * x) * np.exp(-1 * (np.pi**2) * (1/2 + k)**2 * t)
    return result

# Calcular la derivada de la solución respecto a x para los extremos x=0 y x=1 a lo largo del tiempo
flux_x0 = [-derivative_solution_x(0, t, N) for t in t_values]
flux_x1 = [-derivative_solution_x(1, t, N) for t in t_values]

# Graficar el flujo en los extremos x=0 y x=1 a lo largo del tiempo
plt.figure(figsize=(8, 5))
plt.plot(t_values, flux_x0, label='Flujo en x=0', color='blue')  # Cálculo del flujo en x=0
plt.plot(t_values, flux_x1, label='Flujo en x=1', color='red')   # Cálculo del flujo en x=1
plt.xlabel('Tiempo [s]')
plt.ylabel(r'Flujo de Calor [$\frac{^\circ C}{m^2 \cdot s}$]')
plt.title('Flujo de Calor en los Extremos de la Varilla ')
plt.legend()
plt.grid(True)
plt.show()

Observamos que todo el calor aportado en el instante inicial se disipa completamente por la izquierda de la varilla, ya que el flujo en [math] x=1[/math] es nulo, verificando así las condiciones fronteras de nuestro problema.

7 Solución fundamental

La solución fundamental de la ecuación del calor es una función que describe la distribución de temperatura en un medio en función del tiempo y la posición. Para la ecuación del calor unidimensional, la solución fundamental está dada por la expresión:

[math] \phi(x, t) = \frac{1}{2 \sqrt{\pi t}} \exp\left(-\frac{x^2}{t}\right) [/math]

Además la solución general del problema de Cauchy es la solución fundamental (sin tener en cuenta la condición inicial). El problema de Cauchy es el siguiente. Dada la ecuación diferencial parcial:

[math] u_t = k \cdot u_{xx}, [/math]

junto con la siguiente condición inicial:

[math] u(x, 0) = f(x) [/math]

donde \( f(x) \) es una función dada en el dominio \( 0 \leq x \leq L \) y \( k \) es una constante positiva que representa la difusividad térmica del medio.

Para resolver el problema de Cauchy, podemos utilizar la solución fundamental \( \phi(x, t) \) en combinación con la función inicial \( f(x) \) para obtener la solución general \( u(x, t) \) en cualquier instante de tiempo \( t \). Esto se logra convolucionando la función inicial \( f(x) \) con la solución fundamental \( \phi(x, t) \) . La solución \( u(x, t) \) del problema de Cauchy se puede expresar como:

[math] u(x, t) = \int_{-\infty}^{\infty} \phi(x - y, t) f(y) \, dy [/math]

En primer lugar, procederemos a visualizar una representación tridimensional de la solución fundamental \( \phi(x, t) \).

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

# Definir la función que representa la solución fundamental de la ecuación del calor
def phi(x, t):
    return 1 / (2 * np.sqrt(np.pi * t)) * np.exp(-x**2 / t)

# Definir los valores de x y t para la gráfica
x_values = np.linspace(-1, 1, 100)    # Valores de x
t_values = np.linspace(0.01, 1, 200)  # Valores de t (tiempo)

# Crear una malla de valores de x y t
X, T = np.meshgrid(x_values, t_values)

# Calcular la solución fundamental para cada par de valores (x, t)
Z = phi(X, T)

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

# Graficar la superficie 3D
surf = ax.plot_surface(X, T, Z, cmap='viridis')

# Etiquetas y título con el símbolo de "phi"
ax.set_xlabel(r'$x$')
ax.set_ylabel(r'$t$')
ax.set_zlabel(r'$\phi(x, t)$')
ax.set_title(r'Solución fundamental de la ecuación del calor: $\phi(x, t)$')

# Añadir una barra de colores para indicar los valores de Z
fig.colorbar(surf, shrink=0.5, aspect=8)
plt.show()    # Mostrar la gráfica

En esta ocasión veremos la solución de la ecuación del calor \( u_t - u_{xx} = 0 \) en \( x \in \mathbb{R} \) asociada al dato inicial \( f(x) = 1_{[-1,1]} \), viene dada por la convolución:

[math] u(x, t) = \int_{-\infty}^{\infty} \phi(x - y, t) f(y) \, dy = \int_{-\infty}^{\infty} \frac{1}{2 \sqrt{\pi t}} \exp\left(-\frac{(x-y)^2}{t}\right) 1_{[-1,1]} \, dy [/math]

Dibujaremos la solución en diferentes instantes de tiempo: \( t = 0.001 \), \( t = 0.01 \) , \( t = 0.1 \), \( t = 1 \). Aproximaremos la integral indefinida utilizando la fórmula del trapecio.

import numpy as np
import matplotlib.pyplot as plt

# Definir la solución fundamental
def phi(x, t):
    return 1 / np.sqrt(4 * np.pi * t) * np.exp(-x**2 / (4 * t))

def u0(x):
    return np.where(np.logical_and(x >= -1, x <= 1), 1, 0)

# Definir el intervalo espacial y los instantes de tiempo
x = np.linspace(-5, 5, 1000)  # Valores de x
t_values = [0.001, 0.01, 0.1, 1]  # Valores de t (instantes de tiempo)

# Crear una nueva figura para la gráfica
plt.figure(figsize=(10, 6))

# Calcular y graficar la solución para diferentes instantes de tiempo
for t in t_values:
    u_x_t = []
    for xi in x:
        # Utilizar np.trapz para calcular la solución integral de la ecuación del calor
        u_x_t.append(np.trapz(phi(xi - x, t)*u0(x), x))
    plt.plot(x, u_x_t, label=f't = {t}')

# Etiquetas y título de la gráfica
plt.xlabel('x')
plt.ylabel('u(x, t)')
plt.title(f'Solución de $u_t - uxx = 0$ en diferentes instantes de tiempo')
plt.legend()
plt.show()

Por último, la metodología para calcular la solución fundamental de la ecuación del calor es flexible cuando se trata de incorporar condiciones de contorno. Supongamos que deseamos calcular una solución para la ecuación del calor en una dimensión en el semiespacio donde \( x > 0 \), satisfaciendo las condiciones \( u(0, t) = 1 \), \( u(x, t) \rightarrow 0 \) cuando \( x \rightarrow \infty \) y \( u(x, 0) = 0 \) para \( x > 0 \). Podemos modificar la solución fundamental de manera que satisfaga estas condiciones. Al hacerlo, obtenemos la expresión \( u(x,t) = -\text{erf}\left(\frac{x}{2 \cdot \sqrt{t}}\right) + 1 \), donde \( \text{erf}(x) = \frac{2}{\sqrt{\pi}} \int_{0}^{x} e^{-z^2} \, dz \). Esta expresión ajustada permite obtener la solución deseada que cumple con las condiciones de contorno especificadas.

# Importar las bibliotecas necesarias
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.special import erf  # Importar la función de error gaussiano

# Definir la función u(x, t)
def u(x, t):
    return -1 * erf(x / (2 * np.sqrt(t))) + 1

# Definir los valores de x y t para la gráfica
x_values = np.linspace(0, 5, 500)  # Valores de x
t_values = np.linspace(0.01, 5, 500)  # Valores de t (tiempo)

# Crear una malla de valores de x y t
X, T = np.meshgrid(x_values, t_values)

# Calcular los valores de u(x, t) para cada par de valores (x, t)
Z = u(X, T)

# Crear una figura en 3D
fig = plt.figure(figsize=(10, 6))
ax = fig.add_subplot(111, projection='3d')
surf = ax.plot_surface(X, T, Z, cmap='viridis')

ax.set_xlabel('x')
ax.set_ylabel('t')
ax.set_zlabel('u(x, t)')
ax.set_title('Función u(x, t)')
plt.show()

7.1 Solución fundamental dimensión 2

Veamos el problema de Cauchy para la ecuación del calor en dos dimensiones (\( x_1, x_2 \) y \( t \)). Sea \( u(x_1, x_2, t) \) la función de temperatura en un medio bidimensional, que satisface la ecuación del calor:

[math] u_t = k(u_{x_1 x_1} + u_{x_2 x_2}) [/math]

Con las condiciones iniciales dadas por la temperatura inicial \( u(x_1, x_2, 0) = f(x_1, x_2) \) para algún campo de temperatura dado \( f(x_1, x_2) \text{para} (x_1, x_2) \in \partial \Omega \), donde \( \Omega \) es el dominio bidimensional sobre el cual se define la temperatura. La solución general a este problema (sin contar la condición inicial) viene dada por la siguiente expresión:

[math] \phi(x_1, x_2, t) = \frac{1}{4 \pi t} \exp\left( -\frac{x_1^2 + x_2^2}{4t} \right) [/math]

Procederemos ahora a representar la solución fundamental, donde veremos un gráfico tridimensional, [math] x_1, x_2, \phi(x_1, x_2, t) [/math] donde se ha fijado un tiempo t.

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

# Definir la solución fundamental en 2D
def phi_2D(x1, x2, t):
    if t == 0:
        return np.zeros_like(x1)  # Devuelve una matriz de ceros si t es cero
    else:
        return 1 / (4 * np.pi * t) * np.exp(-(x1**2 + x2**2) / (4 * t))

# Definir el intervalo espacial
x1_values = np.linspace(-1, 1, 100)  # Valores de x1
x2_values = np.linspace(-1, 1, 100)  # Valores de x2
x1, x2 = np.meshgrid(x1_values, x2_values)  # Malla de valores de x1 y x2

# Definir los tiempos
t_values = np.linspace(0, 0.1, 3)  # Valores de tiempo

# Tamaño de la figura
fig = plt.figure(figsize=(10, 6))

# Graficar la solución fundamental en 2D para cada tiempo
for i, t in enumerate(t_values, 1):
    ax = fig.add_subplot(1, len(t_values), i, projection='3d')
    phi = phi_2D(x1, x2, t)
    ax.plot_surface(x1, x2, phi, cmap='viridis')
    ax.set_title(f'Tiempo $t=${round(t, 3)}')
    ax.set_xlabel('x1')
    ax.set_ylabel('x2')
    ax.set_zlabel('phi(x1, x2, t)')
plt.show()