Revisión del 11:17 7 mar 2024

1 . Introducción

La visualización de soluciones de la ecuación del calor en una dimensión constituye una herramienta poderosa para desentrañar los intricados patrones de propagación térmica en sistemas unidimensionales. En este estudio, nos sumergimos en la riqueza matemática que rodea la representación gráfica de soluciones, explorando la diversidad de comportamientos térmicos que emanan de este modelo fundamental. A través de la creación visual, buscamos capturar la esencia dinámica de la ecuación del calor y ofrecer una perspectiva única sobre cómo la temperatura evoluciona en función del tiempo y la posición.

La ecuación del calor en una dimensión sirve como un punto de partida teórico sólido para comprender la difusión térmica en estructuras lineales. Al trazar diferentes soluciones, nos sumergimos en el tejido mismo de la propagación térmica, revelando detalles intrincados que, de otra manera, podrían escapar a una descripción puramente analítica. Este enfoque gráfico no solo ilustra la dinámica temporal de la temperatura, sino que también destaca la influencia de condiciones iniciales y parámetros en la evolución térmica a lo largo del espacio unidimensional.

2 . Planteamiento del sistema

Se considera una varilla metálica que ocupa el intervalo [0, 1] y que se encuentra aislada por su superficie lateral, de manera que la conducción de calor sólo se produce en la dirección longitudinal. La temperatura inicial de la varilla es 0 grados. En el extremo izquierdo se consigue mantener la temperatura a 0 grados mientras que en el derecho la temperatura es siempre de 1 grado. De acuerdo con todas estas indicaciones, el sistema de EDP que modeliza el comportamiento de la temperatura es el siguiente:

donde [math]u(x,t)[/math] es la temperatura en la posición [math]x[/math] en el tiempo [math]t[/math], [math]\kappa\gt0[/math] es la constante de conductividad térmica y [math]c[/math] es el calor específico del material, el cual también se considera constante. Las constantes [math]\kappa[/math] y [math]c[/math] provienen de las relaciones:

donde [math]q[/math] es el flujo de calor y [math]e[/math] es la energía calorífica. Cabe destacar que por simplificar el problema se suelen tomar los valores [math]\kappa[/math] y [math]c[/math] como constantes tal y como se ha hecho en este caso, ya que en bastantes casos de interés es razonable ignorar su variación. Si simplificamos aún más, suponiendo que ambas constantes son igual a 1, se tiene el siguiente sistema:

3 . Deducción de la Solución estacionaria

Antes de tratar de resolver el sistema, intentemos conjeturar qué podría suceder. Dado que la temperatura en el lado derecho es mayor que la del lado izquierdo según las condiciones frontera, el calor comenzará a fluir desde el extremo más caliente, aumentando la temperatura dentro de la barra y causando una salida de calor hacia el lado frío. Por otro lado, el aumento interior de la temperatura hace que el flujo de entrada de calor disminuya con el tiempo (ya que el flujo de entrada es exactamente la derivada de la temperatura con respecto a x cambiada de signo), mientras que el flujo de salida aumenta. Esperamos que tarde o temprano los dos flujos se equilibren entre sí y que la temperatura eventualmente alcance una distribución de estado estacionario. Por tanto, es razonable suponer que si se hace el límite cuando el tiempo tiende a infinito de la función temperatura [math]u[/math], se obtendrá una función que tan sólo depende de [math]x[/math], ya que la temperatura debe haber llegado al estado estacionario mencionado en algún momento. Por tanto, bajo esta suposición, la siguiente función [math]v[/math] está bien definida (depende sólo de x):

Además, cumple el sistema de EDO:

Si resolvemos este problema de valor inicial obtendremos:

Es decir, que una vez que la temperatura se estabilice y deje de variar con el tiempo, seguirá el comportamiento de esta función lineal [math]v[/math]. Dicho resultado será muy importante ya que nos permitirá homogeneizar el sistema planteado en la anterior sección.

4 . Resolución de la parte no estacionaria

Para comenzar a resolver el sistema de EDP inicial, homogeneizamos las condiciones frontera utilizando la solución estacionaria que hemos obtenido en el anterior apartado, planteando un problema equivalente con condiciones de frontera homogéneas. Para ello, es conveniente introducir la función:

De esta forma, la solución se compondrá de dos partes, la parte estacionaria y la parte no estacionaria, es decir que la podemos escribir de la siguiente forma:

donde [math]w(x,t)[/math] corresponde a la parte no estacionaria de la solución que es la abordaremos en esta sección. Veamos como podemos obtener las nuevas condiciones frontera y condiciones iniciales:

La función [math]w[/math] introducida nos permite llegar al sistema homogéneo, que será de la forma:

Ahora procedemos a la resolución del sistema por medio del método de separación de variables, es decir suponiendo que la solución es de la forma [math]w(x,t)=X(x)T(t)[/math]. De esta suposición, se obtienen dos problemas de valor inicial uno asociado a la parte espacial de la solución X(x) y el otro asociado a la parte temporal T(t). El primero de ellos será:

Para este problema tendremos como solución cualquiera de las siguientes auto-funciones [math]X_k=B\sin{(k \pi x)} [/math] para [math]k \in N [/math], con [math]B[/math] constante. Por otro lado, para el problema [math] T'+ \lambda T= 0 [/math] obtendremos como solución cualquiera de las siguientes auto-funciones [math] T_k= Ce^{-k^2\pi^2t} [/math] para [math]k \in N [/math], con [math]C[/math] constante. Sin embargo, el producto entre cualquiera de estas auto-funciones nos da una función solución de la ecuación diferencial y las condiciones frontera, que sin embargo, no cumple con la condición inicial. Para tratar de conseguir una solución que cumpla también la condición inicial proponemos la siguiente solución resultante de hacer una suma infinita del producto de las auto-funciones y hacer el desarrollo en serie de Fourier de la condición inicial del sistema para elegir los coeficientes que multiplican a los productos:

De modo que la expresión para la solución final será:

Es fácil ver que la serie converge uniformemente en [math][0,1]\times [s,\infty), \forall s\gt0[/math] mediante el criterio de Weierstrass, al igual que las series de las derivadas. Por lo tanto, es posible derivar la solución metiendo las derivadas dentro de la serie de forma que cumple la ecuación diferencial. Además, como la serie converge uniformemente en los extremos de las condiciones frontera también es posible meter el límite cuando [math]x[/math] tiende tanto a 0 como a 1 dentro de la serie obteniendo así que se cumplen las condiciones frontera. Estas dos comprobaciones son necesarias a la hora de calcular el flujo de calor en los extremos más adelante, por ello se han incluido al final de esta sección. Sin embargo, si se cambia el sistema de EDP, no se repetirá de nuevo esta comprobación y simplemente se dará por hecho (el lector puede comprobarlo si así desea).

5 . Dibujo para [math]t\in[0,1][/math]

En esta sección veremos visualmente cuál es la forma de la solución obtenida anteriormente. Para ello, aproximaremos la forma de la solución considerando los 10 primeros términos de la serie. Además, tan sólo la dibujaremos en el intervalo temporal [math][0,1][/math], ya que este breve intervalo de tiempo es suficientemente grande como para llegar a observar claramente cómo la solución se acerca a su estado estacionario. Representarla para tiempos mayores que este no tiene sentido, teniendo en cuenta que la solución seguirá manteniéndose estacionaria con el tiempo tal y cómo ya sabemos.

Una vez dicho esto, dibujamos la aproximación de la solución u de forma tridimensional:

Aproximación de la función u, dibujada para [math](x,t)\in[0,1]\times[0,1].[/math]

Se puede ver claramente cómo a medida que avanza el tiempo la solución va llegando a un estado estacionario. Esto se ve ya que la solución comienza a dejar de variar con el tiempo. Además, este estado estacionario viene claramente dado por la solución estacionaria obtenida anteriormente [math]v(x)=x[/math]. En la gráfica, se puede ver cómo a medida que pasa el tiempo la solución se va aproximando a dicha solución estacionaria la cual al estar representada tridimensionalmente tiene forma de plano. Además, podemos observar algunas oscilaciones cuando [math]t=0[/math] resultantes de hacer la serie de Fourier en la condición inicial. Esto es debido al fenómeno de Gibbs. La extensión periódica de la extensión impar en [math][-1,1][/math] de la función [math]f(x)=-x[/math] no es continua. En particular, en el punto 1 no lo es, por tanto, a la hora de hacer la serie de Fourier se produce el fenómeno de Gibbs alrededor del 1.

5.1 . Código

Este código realiza una aproximación a la solución de una ecuación del calor unidimensional utilizando la expansión en series de Fourier. Se generan las matrices correspondientes y se grafica la solución en un espacio tridimensional.

% Limpia las variables y la ventana de figuras
clear all
close all

% Define la función de calor inicial f(x) = x
f = @(x) x;

% Genera puntos en el dominio espacial y temporal
xx = 0:0.01:1;
tt = 0:0.01:1;

% Inicializa la matriz A para almacenar los términos de la serie de Fourier
A = zeros(10, length(xx));

% Calcula los términos A(k,:) de la serie de Fourier
for k = 1:10
    A(k, :) = 2 * (-1)^k / (k * pi) * sin(k * pi * xx);
end

% Inicializa la matriz B para almacenar los términos de la serie de Fourier en el tiempo
B = zeros(10, length(tt));

% Calcula los términos B(k,:) de la serie de Fourier en el tiempo
for k = 1:10
    B(k, :) = exp(-k^2 * pi^2 * tt);
end

% Inicializa la matriz Suma para almacenar la suma de los términos de la serie de Fourier
Suma = zeros(length(xx), length(tt));

% Calcula la suma de los términos de la serie de Fourier
for k = 1:10
    Suma = Suma + A(k, :)' * B(k, :);
end

% Calcula la solución u(x,t) sumando la función inicial y la suma de la serie de Fourier
U = repmat(f(xx)', 1, length(tt)) + Suma;

% Gráficas
surf(xx, tt, U', 'FaceColor', 'interp')
xlabel('x')
ylabel('t')
zlabel('u(x,t)')
legend('u(x,t)')

6 . Flujo entrante y saliente

Una vez nos hemos familiarizado con la forma de la solución, veamos cómo se comporta el flujo del calor en los dos extremos de la varilla. Tal y como se dijo anteriormente, el flujo del calor viene dado por:

donde el segundo igual se debe a la elección de [math]\kappa=1[/math]. Anteriormente vimos que era posible derivar la serie tan sólo derivando lo de dentro, de forma que:

Por lo tanto, tendremos que el flujo de calor vendrá dado por:

En este caso, como queremos observar el flujo de calor en los extremos de la varilla, tendremos que el flujo en un tiempo [math]t[/math] en el extremo [math]x=0[/math] vendrá dado por la función de [math]t[/math], [math]q(0,t)[/math] y el flujo en un tiempo [math]t[/math] en el extremo [math]x=1[/math] vendrá dado por la función de [math]t[/math], [math]q(1,t)[/math]. Estas funciones tienen las siguientes expresiones:

Si dibujamos ambas en función del tiempo para [math]t\in[0,1][/math], tenemos lo siguiente:

Aproximación de la función u, dibujada para [math](x,t)\in[0,1]\times[0,1].[/math]

En la primera gráfica podemos ver cómo el flujo de calor comienza siendo nulo y empieza a hacerse cada vez más negativo con el tiempo hasta un punto en el cual se mantiene igual a -1. Como el flujo de calor tiene signo positivo cuando va en la dirección positiva de la [math]x[/math], el que se haga cada vez más negativo significa que cada vez sale más calor por el extremo de la varilla de [math]x=0[/math]. Además, el flujo de calor que sale va aumentando hasta ser igual a -1, punto en el cual se mantiene siempre igual con el tiempo. Esto cuadra con la solución estacionaria obtenida [math]v(x)=x[/math], la cual cumple que [math]-\frac{\partial v}{\partial x}(x)=-1[/math]. Es decir, que el flujo de calor una vez se llega a la solución estacionaria se mantiene constantemente igual a -1 a lo largo de toda la varilla. Por otro lado, en el otro extremo de la varilla, podemos ver como el flujo de calor comienza siendo negativo y aumenta con el tiempo hasta llegar a ser igual a -1 tal y como dicta la solución estacionaria. Es decir, en este caso, empieza entrando calor por este extremo de la varilla (ya que el signo negativo significa que el flujo de calor va en la dirección negativa de las x), aunque el flujo de calor entrante cada vez es menor hasta que se mantiene constantemente igual a -1 cuando se llega a la solución estacionaria.

6.1 .Código

Este código realiza una aproximación a los flujos de calor en los puntos x = 0 y x = 1 utilizando la expansión en series de Fourier. Se calculan los términos correspondientes y se grafican los flujos de calor en función del tiempo.

% Limpia las variables y la ventana de figuras
clear all
close all

% Define la función f_x(x) = 1
f_x = @(x) 1;

% Genera puntos en el tiempo
tt = 0:0.01:1;

% Inicializa la matriz B0 para almacenar los términos de la serie de Fourier en el tiempo para x = 0
B0 = zeros(10, length(tt));

% Calcula los términos B0(k,:) de la serie de Fourier en el tiempo para x = 0
for k = 1:10
    B0(k, :) = 2 * (-1)^k * exp(-k^2 * pi^2 * tt);
end

% Inicializa la matriz B1 para almacenar los términos de la serie de Fourier en el tiempo para x = 1
B1 = zeros(10, length(tt));

% Calcula los términos B1(k,:) de la serie de Fourier en el tiempo para x = 1
for k = 1:10
    B1(k, :) = 2 * (-1)^(2*k) * exp(-k^2 * pi^2 * tt);
end

% Calcula los flujos de calor q(0,t) y q(1,t)
q0 = -f_x(0) - sum(B0);
q1 = -f_x(1) - sum(B1);

% Gráficas
subplot(1, 2, 1)
plot(tt, q0)
xlabel('t')
ylabel('q(0,t)')
legend('q(0,t)')

subplot(1, 2, 2)
plot(tt, q1)
xlabel('t')
ylabel('q(1,t)')
legend('q(1,t)')

7 . Cambio del coeficiente de conductividad térmica

En esta sección, cambiaremos el valor del coeficiente de conductividad térmica a [math]\kappa=1/2[/math] y veremos qué cambios hay en la solución. Repitiendo todo el proceso anterior con este nuevo coeficiente se llega a que la nueva solución es:

Con el fin de dibujar la solución de forma aproximada consideraremos de nuevo los 10 primeros términos de la serie. Además, tan sólo la dibujaremos en el intervalo temporal [math][0,1][/math] de nuevo, ya que este breve intervalo de tiempo es suficientemente grande como para llegar a observar claramente cómo la solución se acerca a su estado estacionario.

Una vez dicho esto, dibujamos la aproximación de la solución u de forma tridimensional:

Aproximación de la función u, dibujada para [math](x,t)\in[0,1]\times[0,1].[/math]

Podemos ver que la forma de la solución es muy parecida a la correspondiente a la elección [math]\kappa=1[/math]. Se ve claramente que se sigue acercando a la misma solución estacionaria a medida que transcurre el tiempo. Además, siguen apareciendo las oscilaciones en la condición inicial provocadas por el fenómeno de Gibbs. Aunque en la gráfica de la solución no se aprecia mucho, en la gráfica de las diferencias entre la solución en el punto medio de la varilla [math]u(1/2,t)[/math] y la solución estacionaria en dicho punto medio [math]v(1/2)[/math] a lo largo del tiempo, podemos observar cómo la nueva solución con [math]\kappa=1/2[/math] tarda más tiempo en llegar al valor de la solución estacionaria en el punto medio de la varilla (ya que la diferencia tarda más tiempo en hacerse 0). Esto no sólo ocurre en el punto medio, sino a lo largo de toda la varilla. La nueva solución tarda más en llegar al estado estacionario (lo cual tiene sentido ya que se tiene un coeficiente de conductividad térmica más pequeño y el flujo de calor entonces será menor en valor absoluto, debido a la fórmula [math]q=-\kappa \frac{\partial u}{\partial x}[/math]). Es decir, cuanto mayor sea el coeficiente de conductividad térmica, menos se tarda en llegar al estado estacionario en la varilla.

7.1 . Código

Este código realiza una aproximación a la solución de una ecuación del calor unidimensional utilizando la expansión en series de Fourier. Se generan las matrices correspondientes y se grafica la solución para dos casos diferentes de [math] \kappa [/math]. Además, se muestra la diferencia entre la solución en [math] x = \frac{1}{2}[/math] la función [math]f(x) [/math] para ambos casos.

% Limpia las variables y la ventana de figuras
clear all
close all

% Define la función f(x) = x
f = @(x) x;

% Genera puntos en el dominio espacial y temporal
xx = 0:0.01:1;
tt = 0:0.01:1;

% Inicializa la matriz A para almacenar los términos de la serie de Fourier
A = zeros(10, length(xx));

% Calcula los términos A(k,:) de la serie de Fourier
for k = 1:10
    A(k, :) = 2 * (-1)^k / (k * pi) * sin(k * pi * xx);
end

% Inicializa las matrices B1 y B2 para almacenar los términos de la serie de Fourier en el tiempo
B1 = zeros(10, length(tt));
B2 = zeros(10, length(tt));

% Calcula los términos B1(k,:) y B2(k,:) de la serie de Fourier en el tiempo
for k = 1:10
    B1(k, :) = exp(-k^2 * pi^2 * tt);
    B2(k, :) = exp(-k^2 * pi^2 / 2 * tt);
end

% Inicializa las matrices Suma1 y Suma2 para almacenar la suma de los términos de la serie de Fourier
Suma1 = zeros(length(xx), length(tt));
Suma2 = zeros(length(xx), length(tt));

% Calcula la suma de los términos de la serie de Fourier para B1
for k = 1:10
    Suma1 = Suma1 + A(k, :)' * B1(k, :);
end

% Calcula la suma de los términos de la serie de Fourier para B2
for k = 1:10
    Suma2 = Suma2 + A(k, :)' * B2(k, :);
end

% Calcula la solución u(x,t) sumando la función inicial y la suma de la serie de Fourier para B1 y B2
U1 = repmat(f(xx)', 1, length(tt)) + Suma1;
U2 = repmat(f(xx)', 1, length(tt)) + Suma2;

% Gráficas
subplot(2, 2, [1, 2])
surf(xx, tt, U1', 'FaceColor', 'interp')
xlabel('x')
ylabel('t')
zlabel('u(x,t)')
legend('u(x,t)')

% Diferencia u(1/2,t) - v(1/2) para B1
subplot(2, 2, 3)
plot(tt, U1(51, :) - f(1/2))
xlabel('t')
ylabel('u(1/2,t) - v(1/2)')
legend('\kappa = 1')

% Diferencia u(1/2,t) - v(1/2) para B2
subplot(2, 2, 4)
plot(tt, U2(51, :) - f(1/2))
xlabel('t')
ylabel('u(1/2,t) - v(1/2)')
legend('\kappa = 1/2')

8 . Cambio de condición inicial y condición frontera

En esta sección supondremos ahora que la temperatura en el lado derecho también es de 0 grados. Además, cambiaremos la condición inicial, de forma que se tendrá el siguiente sistema de EDP que modeliza el comportamiento de la temperatura:

En este caso, como las condiciones frontera ya son homogéneas no es necesario hacer ningún cambio de variable, ya que la solución estacionaria obviamente será la función nula [math]v(x)=0, 0\leq x\leq 1[/math]. Directamente, aplicamos el método de separación de variables llegando a que la solución será:

donde [math]c_{k}[/math] se obtendrán al imponer la condición inicial:

Estos [math]c_{k}[/math] son los coeficientes de la serie de Fourier de la función [math]max \{0,1-4|x-1/2|\}[/math] en [math][0,1][/math] obtenidos tras realizar su extensión impar a [math][-1,1][/math]. Estos coeficientes no se muestran analíticamente ya que se calcularán usando la fórmula del trapecio a la hora de dibujar la función solución [math]u[/math]. Vendrán dados por:

A continuación, dibujaremos la gráfica de la solución. Para ello, de nuevo, aproximaremos la forma de la solución considerando los 10 primeros términos de la serie. Además, tan sólo la dibujaremos en el intervalo temporal [math][0,1][/math], ya que este breve intervalo de tiempo es suficientemente grande como para llegar a observar claramente cómo la solución se acerca a su estado estacionario, al igual que en el caso anterior.

Una vez dicho esto, dibujamos la aproximación de la solución u de forma tridimensional:

Cambiar foto</math>

Se puede ver claramente cómo a medida que avanza el tiempo la solución va llegando a un estado estacionario. Esto se ve ya que la solución comienza a dejar de variar con el tiempo. Además, este estado estacionario viene claramente dado por la solución estacionaria obtenida anteriormente de forma trivial [math]v(x)=0[/math]. En la gráfica, se puede ver cómo a medida que pasa el tiempo la solución se va aproximando a dicha solución estacionaria. Además, podemos observar menos oscilaciones cuando [math]t=0[/math] que en el caso anterior, de hecho si se aumenta el número de términos de la serie para aproximar la solución las oscilaciones prácticamente desaparecen. Esto se debe a que cuando se impone la condición inicial en este caso, la extensión periódica de la extensión impar en [math][-1,1][/math]de la función [math]max \{0,1-4|x-1/2|\}[/math] es continua, de forma que no se produce el fenómeno de Gibbs como ocurría en el caso anterior.

Por otro lado, podemos ver como en [math]x=0[/math] el flujo de calor empieza siendo negativo y se vuelve más negativo con el tiempo hasta cierto punto, a partir del cual empieza a crecer hasta hacerse cero. Es decir, como la dirección positiva del flujo de calor es la dirección positiva de la x, esto significa que el flujo de calor comienza saliendo del lado izquierdo de la varilla cada vez más hasta un punto a partir del cual cada vez sale menos calor hasta que deja de salir cuando la temperatura de la varilla se estabiliza. De igual forma, cuando [math]x=1[/math] el flujo de calor empieza siendo positivo y se vuelve más positivo con el tiempo hasta cierto punto, a partir del cual empieza a crecer hasta hacerse cero. Es decir, el flujo de calor comienza saliendo del lado derecho de la varilla cada vez más hasta un punto a partir del cual cada vez sale menos calor hasta que deja de salir cuando la temperatura de la varilla se estabiliza. En ambos extremos ocurre lo mismo. Esto se debe a que inicialmente en tiempo 0, la condición inicial dada [math]u(x,0)=max \{0,1-4|x-1/2|\}[/math] impone que haya una temperatura mayor que 0 en la parte central de la varilla. Esto hace que llegue una 'ola' de calor hacia ambos lados la cual hace que el flujo de calor saliente por ambos lados aumente. Una vez dicha transmisión de calor se va estabilizando, el flujo de calor saliente va disminuyendo hasta que deja de salir calor, estabilizándose la temperatura siendo igual a 0 grados (la solución estacionaria) a lo largo de toda la varilla.

8.1 . Código

Este código realiza una aproximación a la solución de una ecuación del calor unidimensional utilizando la expansión en series de Fourier. Se calculan los coeficientes de la serie de Fourier, se generan las matrices correspondientes y se grafica la solución y los flujos de calor en los extremos de la barra.

% Limpia las variables y la ventana de figuras
clear all
close all

% Define la función de calor inicial f(x)
f = @(x) max(0, 1 - 4 * abs(x - 1/2));

% Genera puntos en el dominio espacial y temporal
xx = 0:0.01:1;
tt = 0:0.01:1;

% Inicializa el vector a para almacenar los coeficientes de la serie de Fourier
a = zeros(1, length(xx));

% Calcula los coeficientes a(k) utilizando la regla del trapecio
for k = 1:10
    a(k) = 2 * trapz(xx, f(xx) .* sin(k * pi * xx));
end

% Inicializa las matrices A y B para almacenar los términos de la serie de Fourier
A = zeros(10, length(xx));
B = zeros(10, length(tt));

% Calcula las matrices A y B
for k = 1:10
    A(k, :) = a(k) * sin(k * pi * xx);
    B(k, :) = exp(-k^2 * pi^2 * tt);
end

% Inicializa la matriz U para almacenar la solución aproximada
U = zeros(length(xx), length(tt));

% Calcula la solución aproximada u(x, t) sumando los términos de la serie de Fourier
for k = 1:10
    U = U + A(k, :)' * B(k, :);
end

% Gráficas
subplot(2, 2, [1, 2])
surf(xx, tt, U', 'FaceColor', 'interp')
xlabel('x')
ylabel('t')
zlabel('u(x,t)')
legend('u(x,t)')

% Flujos de calor en los puntos x = 0 y x = 1
tt = 0:0.01:1;

% Calcula los flujos de calor q(0,t) y q(1,t)
B0 = zeros(10, length(tt));
B1 = zeros(10, length(tt));

for k = 1:10
    B0(k, :) = a(k) * k * pi * exp(-k^2 * pi^2 * tt);
    B1(k, :) = a(k) * k * pi * (-1)^k * exp(-k^2 * pi^2 * tt);
end

% Suma los flujos de calor para obtener q(0,t) y q(1,t)
q0 = -sum(B0);
q1 = -sum(B1);

% Gráficas
subplot(2, 2, 3)
plot(tt, q0)
xlabel('t')
ylabel('q(0,t)')
legend('q(0,t)')

subplot(2, 2, 4)
plot(tt, q1)
xlabel('t')
ylabel('q(1,t)')
legend('q(1,t)')

9 . Cambio de condición inicial

En esta sección mantendremos el sistema de la sección anterior pero en este caso cambiando la suposición de que en el lado derecho de la varilla la temperatura es de 0 grados por la de que el lado derecho se encuentra aislado térmicamente. Por lo tanto, cambiaremos las condiciones frontera tipo Dirichlet, de forma que se tendrá el siguiente sistema de EDP que modeliza el comportamiento de la temperatura:

Este cambio, se debe a que cuando se dice que la barra está aislada térmicamente en el extremo derecho, significa que no hay intercambio de calor con el entorno en ese extremo específico. Como se vio anteriormente, el flujo de calor viene dado por [math]q=-\kappa \frac{\partial u}{\partial x}=-\frac{\partial u}{\partial x}[/math], donde la segunda igualdad se debe a que en este caso [math]\kappa=1[/math]. Si igualamos el flujo de calor a 0 en el extremo derecho tenemos que:

Una vez se comprende el planteamiento del problema, pasaremos a su resolución. En este caso, como las condiciones frontera ya son homogéneas no es necesario hacer ningún cambio de variable, ya que la solución estacionaria obviamente será la función nula [math]v(x)=0, 0\leq x\leq 1[/math]. Luego podemos aplicar directamente el método de separación de variables llegando a que la solución será:

Para calcular los coeficientes [math]c_{k}[/math] de forma similar a como se hizo anteriormente, vamos a trabajar con la base [math]\{ \sin{(\frac{\pi}{2} \cdot (2k-1) \cdot x)} \}_{k=1}^{\infty} [/math] de [math]L^2(0,1)[/math]. Es fácil verificar que dicha base es ortogonal:

donde [math]\lt.,.\gt_{L^2(0,1)}[/math] es como se denota el producto escalar en [math]L^2(0,1)[/math].

Por lo tanto, para calcular los coeficientes [math]c_{k}[/math] procedemos de la misma manera, debemos imponer que se cumpla la condición inicial, es decir:

Como la base de [math]L^2(0,1)[/math] mostrada anteriormente es ortogonal, podemos calcular los [math] c_{k} [/math] correspondientes a la solución, que son los coeficientes de Fourier en esta nueva base, de la siguiente forma:

donde el segundo igual viene de que:

Estos coeficientes, al igual que en la anterior sección se calcularán con la fórmula del trapecio a la hora de dibujar la solución.

En esta sección, aunque se dijera anteriormente que se dejarían las comprobaciones de que esta es una solución válida para el lector, se hará dicha comprobación, con el fin de justificar la aplicación del teorema del máximo en la siguiente sección. Para ello, vemos primero que la serie es convergente uniformemente en [math][0,1]\times [s,\infty), \forall s\gt0[/math], acotando los términos de dentro de la siguiente forma:

donde [math](x,t)\in [0,1]\times [s,\infty)[/math] y se ha usado que la función [math]f[/math] está acotada por 1. Por el criterio de Weierstrass, podemos decir que en efecto converge uniformemente. Además, para la serie de las derivadas de lo de dentro se tiene de forma similar que también converge uniformemente. Por lo tanto, esto permite derivar la solución introduciendo las derivadas dentro, lo cual confirma que cumple la ecuación diferencial y además permite meter el límite cuando [math]x[/math] tiende tanto a 0 como a 1 dentro de la serie (en el caso del extremo derecho sería dentro de la serie una vez derivada) confirmando así que se cumplen las dos condiciones frontera. Además, por ser la extensión periódica de [math]f[/math] continua, la serie de Fourier converge puntualmente lo cual garantiza que además la solución cumple la condición inicial y es continua en [math]t=0[/math] también, al unir bien la solución cuando [math]t[/math] tiende a 0 en ambos extremos cuando [math]x=0[/math] o [math]x=1[/math].

A continuación, dibujaremos la gráfica de la solución. Para ello, de nuevo, aproximaremos la forma de la solución considerando los 10 primeros términos de la serie. Además, tan sólo la dibujaremos en el intervalo temporal [math][0,1][/math], ya que este breve intervalo de tiempo es suficientemente grande como para llegar a observar claramente cómo la solución se acerca a su estado estacionario.

Una vez dicho esto, dibujamos la aproximación de la solución u de forma tridimensional:

Cambiar foto</math>

De nuevo, podemos ver claramente cómo a medida que avanza el tiempo la solución va llegando a un estado estacionario. Además, este estado estacionario viene claramente dado por la solución estacionaria obtenida anteriormente de forma trivial [math]v(x)=0[/math]. Además, como la condición inicial sigue siendo la misma que en la anterior sección, tampoco se produce el fenómeno de Gibbs por la misma razón, la función [math]max \{0,1-4|x-1/2|\}[/math] definida en [math][0,1][/math] se mantiene continua al extenderla periódicamente.

Por otro lado, podemos ver como en [math]x=0[/math] el flujo de calor empieza siendo negativo y se vuelve más negativo con el tiempo hasta cierto punto, a partir del cual empieza a crecer hasta hacerse cero. Es decir, como la dirección positiva del flujo de calor es la dirección positiva de la x, esto significa que el flujo de calor comienza saliendo del lado izquierdo de la varilla cada vez más hasta un punto a partir del cual cada vez sale menos calor hasta que deja de salir cuando la temperatura de la varilla se estabiliza. Es decir, ocurre exactamente lo mismo que ocurría en el problema de la anterior sección. Sin embargo, cuando [math]x=1[/math] el flujo de calor se mantiene igual a 0 para todo tiempo, tal y como dicta la condición frontera del problema. Aun así, en la gráfica de la solución se puede ver que en el extremo derecho la temperatura llega a crecer un poco con el tiempo, hasta que llega a un pico y vuelve a bajar hasta llegar a su estado estacionario. Esto se debe a que inicialmente en tiempo 0, la condición inicial dada [math]u(x,0)=max \{0,1-4|x-1/2|\}[/math] impone que haya una temperatura mayor que 0 en la parte central de la varilla. Esto hace que llegue la 'ola' de calor hacia ambos lados que se producía en el anterior problema la cual hace que el flujo de calor saliente por el lado izquierdo aumente, debido a que no puede aumentar la temperatura la cual se mantiene en 0. En el otro extremo, debido a que está aislado y no puede haber un flujo de calor saliente, la temperatura aumenta. Una vez dicha transmisión de calor se va estabilizando, el flujo de calor saliente por el extremo izquierdo va disminuyendo hasta que deja de salir calor, mientras que en el extremo derecho la temperatura comienza a bajar, estabilizándose la temperatura de toda la varilla hasta llegar a ser de 0 grados.

9.1 . Código

Este código en MATLAB realiza una aproximación a la solución de una ecuación de onda unidimensional mediante la expansión en series de Fourier. A continuación, te proporciono comentarios detallados para cada sección del código:

% Limpia las variables y la ventana de figuras
clear all
close all

% Define la función de la condición inicial f(x)
f = @(x) max(0, 1 - 4 * abs(x - 1/2));

% Genera puntos en el dominio espacial y temporal
xx = 0:0.01:1;
tt = 0:0.01:1;

% Inicializa el vector c para almacenar los coeficientes de la serie de Fourier
c = zeros(1, length(xx));

% Calcula los coeficientes c(k) utilizando la regla del trapecio
for k = 1:10
    c(k) = 2 * trapz(xx, f(xx) .* sin(pi/2 * (2*k - 1) * xx));
end

% Inicializa las matrices A y B para almacenar los términos de la serie de Fourier
A = zeros(10, length(xx));
B = zeros(10, length(tt));

% Calcula las matrices A y B
for k = 1:10
    A(k, :) = c(k) * sin(pi/2 * (2*k - 1) * xx);
    B(k, :) = exp(-pi^2/4 * (2*k - 1)^2 * tt);
end

% Inicializa la matriz U para almacenar la solución aproximada
U = zeros(length(xx), length(tt));

% Calcula la solución aproximada u(x, t) sumando los 10 primeros términos de la serie de Fourier
for k = 1:10
    U = U + A(k, :)' * B(k, :);
end

% Gráficas
subplot(2,2,[1,2])
surf(xx, tt, U', 'FaceColor', 'interp')
xlabel('x')
ylabel('t')
zlabel('u(x,t)')
legend('u(x,t)')

% Flujos de calor en los puntos x = 0 y x = 1

% Calcula los términos dentro de la suma implicada en el cálculo de los flujos
B0 = zeros(10, length(tt));
B1 = zeros(10, length(tt));

for k = 1:10
    B0(k, :) = c(k) * pi/2 * (2*k - 1) * exp(-pi^2/4 * (2*k - 1)^2 * tt);
    B1(k, :) = c(k) * pi/2 * (2*k - 1) * 0 * exp(-pi^2/4 * (2*k - 1)^2 * tt);
end

% Suma los términos y añade el signo negativo para obtener q(0,t) y q(1,t)
q0 = -sum(B0);
q1 = -sum(B1);

% Gráficas
subplot(2,2,3)
plot(tt, q0)
xlabel('t')
ylabel('q(0,t)')
legend('q(0,t)')

subplot(2,2,4)
plot(tt, q1)
xlabel('t')
ylabel('q(1,t)')
legend('q(1,t)')

10 . Principio del Máximo

El teorema del principio del máximo nos permite asegurar la unicidad del sistema. Dicho teorema nos dice lo siguiente:

Sea [math] u \in C^{2,1} (Q_T) \cap C(\overline{Q_T})[/math] con [math] Q_T=\Omega \times (0,T)[/math] siendo [math]\Omega[/math] abierto de [math] \mathbb R^n[/math], que verifica:

Entonces [math]u[/math] alcanza su máximo en la frontera parabólica, es decir:

En el caso de nuestra solución correspondiente al apartado anterior, se usaría [math]n=1[/math] y [math]\Omega=(0,1)[/math]. Como la serie converge uniformemente en [math][0,1]\times [s,\infty), \forall s\geq 0[/math] y la serie de las derivadas con respecto a ambas variables también lo hacen (lo vimos en la anterior sección), es posible derivar la función simplemente haciendo la suma de las derivadas de lo de dentro en este dominio tal y como se hizo para obtener el flujo de calor en los extremos. Por lo tanto, como las funciones de dentro son [math]C^{\infty,\infty}((0,1)\times (0,\infty))[/math], la función [math]u[/math] también lo es, debido a que se puede aplicar esto a cada serie de las derivadas. En particular, es [math]C^{\infty,\infty}(Q_T)[/math] para cualquier [math]T\gt0[/math]. Además, como la serie converge uniformemente en las condiciones frontera, es posible meter el límite dentro de la serie, tal y como se mencionó en la anterior sección, de forma que se cumple la continuidad en [math][0,1]\times(0,\infty)[/math]. Faltaría tan sólo ver que la solución [math]u[/math] es continua también para [math][0,1]\times\{0\}[/math]. Esto es fácil de comprobar, ya que como la función de la condición inicial al extenderla periódicamente sigue siendo continua, la serie de Fourier converge puntualmente, además la solución une bien cuando [math]t[/math] tiende a 0 en ambos extremos cuando [math]x=0[/math] o [math]x=1[/math]. Por lo tanto, se concluye que [math]u\in C([0,1]\times[0,T])=C(\overline{Q_T})[/math] para cualquier [math]T\gt0[/math]. Además, por ser solución del sistema de EDPs planteado anteriormente, cumple la desigualdad del principio del máximo al cumplir la igualdad.

Con todo esto, tenemos que cumple todas las hipótesis del teorema del máximo, en primer lugar [math] u \in C^{2,1} (Q_T) \cap C(\overline{Q_T})[/math] para cualquier [math]T\gt0[/math], además por el apartado anterior sabemos que la u que es la solución correspondiente al sistema cumple la expresión [math]\frac{\partial u}{\partial t}- \frac{\partial ^2 u}{\partial x^2}\leq 0 [/math].

10.1 . Unicidad

Sabiendo esto, podemos demostrar la unicidad de la solución haciendo uso del principio del máximo. Si consideramos dos soluciones [math]u[/math] y [math]v[/math] del sistema, si definimos [math]w=u-v[/math] que satisface:

Aplicando el principio del máximo sobre [math]w[/math] obtenemos: [math] \min \limits_{(x,t) \in \partial _P Q_T}w \leq w \leq \max \limits_{(x,t) \in \partial _P Q_T}w [/math]. De esto deducimos que [math]w \equiv 0[/math], y por tanto podemos concluir que [math]u=v[/math], quedando mostrado que la solución del sistema es única.

11 . Solución fundamental de la ecuación del calor en dimensión 1

En este apartado mostraremos cual es el comportamiento de la solución fundamental ecuación del calor de forma visual. Para poder hacer una buena representación, la dibujaremos para [math] x \in [−1, 1] [/math] y [math] t \in [10^{-2}, 1] [/math]. Se define la solución fundamental en una dimensión como:

Su representación es la siguiente:

Cambiar foto</math>

11.1 .Código

clear all
close all

% Define la función solución fundamental de la ecuación del calor
u = @(x, t) exp(-x.^2./(4*t))./sqrt(4*pi*t);

% Define los rangos de x y t
xx=-1:0.01:1;        % Rango de x
tt=10^(-2):0.01:1;   % Rango de t

% Evalúa la función en la malla de puntos (x, t)
[X, T] = meshgrid(xx, tt);
U = u(X, T);

% Grafica la función
surf(X, T, U);
xlabel('x');
ylabel('t');
zlabel('u(x,t)');
title('Solución fundamental de la ecuación del calor');

12 . Solución de la ecuación del calor para el semi-espacio positivo

Continuando con la idea de representar la solución de la ecuación del calor se impondrá una serie de condiciones sobre el sistema para ver como se comporta la solución de acuerdo con el ejercicio 10b) de la hoja 2.

La solución de [math]u_{t} −u_{xx} = 0[/math] en [math] x\gt0, t\gt0[/math] debe satisfacer las condiciones; la priemra es que u(0,t)=1, la segunda es que [math]u(x,t) \rightarrow 0 [/math] si [math]x \rightarrow \infty[/math] para [math] x\gt0 [/math] y la tercera es [math] u(x,0)=0[/math] para [math] x\gt0 [/math]. Vamos a mostrar todas estas condiciones de forma que le sea más sencillo al lector entender que tipo de condiciones se están imponiendo en nuestro sistema.

Para poder realizar este apartado necesitaremos una solución para el problema planteado, por ello usaremos la indicación que se da en el propio ejercicio donde se dice que la solución será de la forma [math]u(x,t)=U(\frac{x}{\sqrt{t}})[/math], que es una solución autosemejante. Por ello antes de resolver el sistema se va a explicar que es una solución autosemejante, la idea es poder entender mejor aún su comportamiento.

Una solución auto-similar o autosemejante se refiere a un tipo particular de solución que exhibe simetría bajo ciertas transformaciones de escala. Esto significa que si ampliamos o reducimos la solución por un factor constante, la forma general de la solución se mantiene inalterada. Este comportamiento auto-similar es especialmente interesante y útil en el estudio de fenómenos que exhiben patrones repetitivos a diferentes escalas. Ahora se procede a la solución del sistema, sustituyendo [math]u(x,t)=U(\frac{x}{\sqrt{t}})[/math] en la ecuación del calor obtenemos lo siguiente:

​ Se trata de una EDO no lineal de orden 2, para poder resolverla debemos realizar un cambio de variable que es [math]y= \frac{x}{\sqrt{t}}[/math] de manera que llegamos a la siguiente expresión:

Se trata de una EDO exacta que si la resolvemos obtendremos la siguiente expresión:

si imponemos las condiciones iniciales especificadas inicialmente se llega a la solución final es:

esta será la expresión que será objeto de nuestro estudio donde se representará y se detallarán los elementos más importantes que la caracterizan.

12.1 .Código

13 . Ecuación del calor y la convolución

Variando las condiciones iniciales de la ecuación del calor podemos obtener diferentes soluciones, en este caso si imponemos que el dato inicial sea [math]u_{0}(x) = 1_[−1,1] [/math] su solución estará definida por la convolución: [math]u(x,t) = \int_{-\infty}^{\infty} K(x-y) \cdot u_{0}(y) dy [/math]

donde K(x,t) es la solución fundamental presentada anteriormente. Se mostrará la solución para diferentes instantes de tiempo, para aproximar la integral que viene dada por la convolución usaremos el método numérico del trapecio.

Cambiar foto</math>

13.1 .Código

clear all
close all

% Define la función solución fundamental de la ecuación del calor
K = @(x,t) exp(-x.^2./(4*t))./sqrt(4*pi*t);

% Define el rango de x,y,t
xx=-5:0.1:5;                % Rango de x
yy=-1:0.01:1;               % Rango de y
tt=[0.001,0.01,0.1];        % Rango de y

%Define la función u
U=zeros(length(tt),length(xx));
for j=1:length(tt)
    for i=1:length(xx)
        U(j,i)= trapz( yy, K( xx(i)*ones(1,length(yy)) - yy, tt(j) ) );
    end
end

% Grafica la función para t = 0.001, t = 0.01 y t = 0.1
subplot(2,3,1)
plot(xx,U(1,:))
xlabel('x');
ylabel('u(x,0.001)');
title('Solución para t=0.001');

subplot(2,3,2)
plot(xx,U(2,:))
xlabel('x');
ylabel('u(x,0.01)');
title('Solución para t=0.01');

subplot(2,3,3)
plot(xx,U(3,:))
xlabel('x');
ylabel('u(x,0.1)');
title('Solución para t=0.1');

%Gráfica de la solución en 3 dimensiones

%Nuevo vector de tiempos
tt=0:0.01:1;
%Define la función u
U=zeros(length(tt),length(xx));
for j=1:length(tt)
    for i=1:length(xx)
        U(j,i)= trapz( yy, K( xx(i)*ones(1,length(yy)) - yy, tt(j) ) );
    end
end

%Gráfica
subplot(2,3,[4,5,6])
surf(xx, tt, U);
xlabel('x');
ylabel('t');
zlabel('u(x,t)');
title('Solución');

14 . Ecuacion del calor de dimensión 2

Ahora se mostrará cual es la ecuación del calor para el caso bidimensional donde se necesitan 3 variables para abordar este sistema, dicha expresión tendrá la forma:

[math]{\partial u\over \partial t} = \alpha \left({\partial^2 u\over \partial x^2 } + {\partial^2 u\over \partial y^2 } + {\partial^2 u\over \partial z^2 }\right)[/math] [math] = \alpha ( u_{xx} + u_{yy} + u_{zz} ) \quad [/math]

donde:

• [math]u = u \left(x, y, z,t\right)[/math] es la temperatura como una función del espacio y del tiempo;
• [math]\frac{\partial u}{\partial t}[/math] es la razón de cambio de temperatura en un punto respecto del tiempo;
• [math]u_{xx},u_{yy}, u_{zz}[/math] son derivadas segundas parciales (conducciones térmica) de la temperatura respecto de [math]x,y,z[/math], respectivamente;
• [math]\alpha = \frac{k}{c_p\rho}[/math] es la difusividad térmica, una cantidad específica del material que depende de la conductividad térmica [math]k[/math], la densidad de masa [math]\rho[/math], y la calor específico [math]c_p[/math].

La ecuación del calor es una consecuencia de la ley de Fourier de conducción . Las soluciones de la ecuación del calor se caracterizan por una suavidad gradual de la distribución de temperatura inicial por el flujo de calor desde las áreas más cálidas hacia las más frías de un objeto. Generalmente, muchos estados diferentes y condiciones iniciales tienden al mismo equilibrio estable. Veremos además que la solución según tiende a t=0 se vuelve mas singular por lo que podemos dar una serie de conclusiones a partir de la representación del sistema.

14.1 .Código

clear all
close all

% Define la función solución fundamental de la ecuación del calor
u = @(x, y, t) exp(-(x.^2+y.^2)./(4*t))./4*pi*t;

% Define los rangos de x, y, t
xx=-1:0.01:1;          % Rango de x
yy=-1:0.01:1;          % Rango de y
tt=[0.001,0.01,0.1];   % Rango de t

% Evalúa la función en la malla de puntos (x, t)
[X, Y] = meshgrid(xx, yy);

for i=1:length(tt)
    U = u(X, Y, tt(i));
    figure(i)
    surf(X, Y, U);
    xlabel('x');
    ylabel('y');
    zlabel("u(x,y,"+num2str(tt(i))+")");
    title("Solución fundamental de la ecuación del calor con t=" + num2str(tt(i)) + ".");
end

Cambiar foto</math>

Cambiar foto</math>

Cambiar foto</math>