Revisión actual del 23:55 7 mar 2024

1 . Introducción

En este trabajo se busca entender con mayor profundidad la ecuación del calor que modeliza la temperatura de una varilla de longitud 1. Además, en la parte final también se estudia la ecuación del calor de una dimensión en [math]\mathbb R[/math] y brevemente la ecuación del calor de dimensión 2 en [math]\mathbb R[/math].

En la primera parte, que es la que ataña a la ecuación del calor que modeliza la varilla, se harán cambios en las condiciones iniciales con el fin de estudiar los diferentes comportamientos que tiene la solución dependiendo de la condición inicial que se ha impuesto. Los cambios en la solución, se relacionarán con el caso real de la varilla y se explicará el sentido físico de estos cambios, observando qué ocurre con los flujos de calor en determinados puntos de la varilla y cómo se estabiliza la temperatura a lo largo del tiempo. También se cambiará el coeficiente de conductividad térmica con el fin de observar si la solución cambia y de qué forma lo hace, relacionando esto con los cambios de temperatura en la varilla. Finalmente, se hará una breve sección teórica hablando sobre el principio del máximo y cómo este es posible aplicarlo en algunos casos para demostrar la unicidad de la solución de un sistema de EDP.

Por otro lado, en la segunda parte del trabajo se hará especial énfasis en la solución fundamental de la ecuación del calor, tanto para dimensión 1 como para dimensión 2. La solución fundamental para dimensión 1 se usará también para hallar una solución que además cumpla una condición inicial específica. Además, en la sección 12, se impondrán algunas condiciones iniciales para las cuales se encontrará una solución 'autosemejante' (más adelante se explicará qué significa esto) que las cumpla todas.

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 solución estacionaria
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 dependientes de x
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 dependientes de t
B = zeros(10, length(tt));

% Calcula los términos B(k,:) de la serie de Fourier
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 solución estacionaria 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:

Gráficas del flujo de calor en [math]x=0[/math] y en [math]x=1[/math]. En la primera se muestra el flujo de calor en [math]\{0\}\times[0,1][/math] en función del tiempo, mientras que en la segunda se muestra el flujo en [math]\{1\}\times[0,1][/math], también en función del tiempo.

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 derivada de la solución estacionaria
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:

En la figura superior se muestra la aproximación de la función u, dibujada para [math](x,t)\in[0,1]\times[0,1][/math]. En las dos figuras inferiores se muestra la diferencia entre la aproximación de la solución en [math]x=1/2[/math] y la solución estacionaria en [math]x=1/2[/math] a lo largo del tiempo para [math]t\in[0,1][/math]. En el caso de la figura inferior izquierda se hace esto para [math]\kappa=1[/math], mientras que en el caso de la figura inferior derecha se hace para [math]\kappa=1/2[/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. Además, grafica las diferencias entre la solución y la solución estacionaria para dos casos diferentes de [math] \kappa [/math] en [math]x=1/2[/math].

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

% Define la función estacionaria
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 dependientes de x
A = zeros(10, length(xx));

% Calcula los términos A(k,:) de la serie de Fourier dependientes de x
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 dependientes de t para cada kappa
B1 = zeros(10, length(tt));
B2 = zeros(10, length(tt));

% Calcula los términos B1(k,:) y B2(k,:) de la serie de Fourier dependientes de t
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 solución estacionaria 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:

En la figura superior se muestra la aproximación de la función u, dibujada para [math](x,t)\in[0,1]\times[0,1][/math]. En las dos figuras inferiores se muestran las gráficas del flujo de calor en [math]x=0[/math] y en [math]x=1[/math]. En la primera se muestra el flujo de calor en [math]\{0\}\times[0,1][/math] en función del tiempo, mientras que en la segunda se muestra el flujo en [math]\{1\}\times[0,1][/math], también en función del tiempo.

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 términos de las sumas de 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 términos y les 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)')

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:

En la figura superior se muestra la aproximación de la función u, dibujada para [math](x,t)\in[0,1]\times[0,1][/math]. En las dos figuras inferiores se muestran las gráficas del flujo de calor en [math]x=0[/math] y en [math]x=1[/math]. En la primera se muestra el flujo de calor en [math]\{0\}\times[0,1][/math] en función del tiempo, mientras que en la segunda se muestra el flujo en [math]\{1\}\times[0,1][/math], también en función del tiempo.

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

% 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\gt 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

A partir de ahora, trataremos de estudiar la ecuación del calor en [math] \mathbb R [/math] en dimensión 1. Más concretamente estudiaremos la EDP con condición inicial:

donde [math]u_0[/math] es una función cualquiera. En esta sección comenzaremos mostrando cual es el comportamiento de la solución fundamental de la ecuación del calor, la cual nos permite obtener soluciones que cumplan la condición inicial deseada. Este comportamiento se mostrará de forma visual, mediante un dibujo de la solución fundamental. 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 recuerda que la solución fundamental en una dimensión se define como:

Su representación es la siguiente:

Representación de la función fundamental u, dibujada para [math](x,t)\in[-1,1]\times[10^{-2},1].[/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 ecuación del calor de dimensión 1 en [math] \mathbb R[/math] se impondrán una serie de condiciones sobre el sistema para ver como se comporta la solución. En este caso, trataremos de encontrar una solución [math]u[/math] tal que:

Para resolver el problema planteado, trataremos de buscar soluciones de la forma [math]u(x,t)=U(\frac{x}{\sqrt{t}})[/math] autosemejantes. Por ello antes de resolver el sistema se va a explicar qué es una solución autosemejante.

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 el cambio de variable [math]y= \frac{x}{\sqrt{t}}[/math], de manera que llegamos a la siguiente expresión:

Esta se trata de una EDO exacta que si se resuelve se obtiene:

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

A continuación, representaremos esta solución en una gráfica. En este caso, se hará la representación para [math](x,t)\in[0,5]\times [0,5][/math], logrando así la siguiente figura:

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

12.1 .Código

clear all
close all

% Define el integrando de la función error
errf = @(z) 2/sqrt(pi)*exp(-z.^2);

% Define los rangos de x y t
xx=linspace(0,5,100);   % Rango de x
tt=linspace(0,5,100);   % Rango de t

%Define la función u
U=zeros(length(tt),length(xx));
for j=1:length(tt)
    for i=1:length(xx)
        limint=linspace(0, xx(i)/( 2*sqrt(tt(j)) ), 100); %Define el intervalo de integración para cada (x,t)
        U(j,i)= 1-trapz( limint , errf(limint) );
    end
end

% Grafica la función
surf(xx, tt, U);
xlabel('x');
ylabel('t');
zlabel('u(x,t)');
title('Solución');

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

Variando las condiciones iniciales de la ecuación del calor podemos obtener diferentes soluciones a partir de la solución inicial. En esta sección, impondremos que el dato inicial sea [math]u_{0}(x) = 1_{[−1,1]}(x) [/math], teniendo así que resolver la EDP con condición inicial:

su solución vendrá dada por la convolución:

donde [math]K[/math] es la solución fundamental presentada anteriormente. Se mostrará la solución para diferentes instantes de tiempo, aproximando la integral que viene dada por la convolución usando el método numérico del trapecio. Haciendo esto, se consigue la siguiente representación:

En la parte superior se muestran representaciones de la solución [math]u[/math] para distintos valores de [math]t[/math] en función de [math]x[/math] para [math]x\in[-5,5][/math]. La primera lo hace para [math]t=0.001[/math], la segunda para [math]t=0.01[/math] y la tercera para [math]t=0.1[/math]. En la parte inferior se muestra la gráfica en 3 dimensiones de la solución para [math](x,t)\in[-5,5]\times[0,1][/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 t

%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 . Solución fundamental de la ecuación del calor en dimensión 2

Para acabar, trataremos de observar la forma de la solución de la ecuación del calor para una dimensión más, dimensión 2. La ecuación del calor para dimensión [math]n\gt1[/math] en [math]\mathbb R[/math] es:

Por lo tanto, la solución fundamental en el caso [math]n=2[/math] es:

Esta función, no es posible representarla de la misma forma que se hizo con la solución fundamental para [math]n=1[/math], ya que al haber dos variables espaciales y una temporal y tener de imagen un escalar, se necesitarían 4 dimensiones para representarla. Sin embargo, podemos fijar varios valores de [math]t[/math] tal y como hicimos anteriormente y representar la gráfica en función de las variables espaciales en 3 dimensiones. Procediendo de esta manera se tiene lo siguiente:

Representación de la solución [math]u[/math] para [math]t=0.1[/math] en función de [math]x[/math] para [math]x\in[-5,5][/math]

Representación de la solución [math]u[/math] para [math]t=0.01[/math] en función de [math]x[/math] para [math]x\in[-5,5][/math]

Representación de la solución [math]u[/math] para [math]t=0.001[/math] en función de [math]x[/math] para [math]x\in[-5,5][/math]

Es posible observar como cuanto más nos acercamos al [math]t=0[/math] más singular es la solución, debido a que la [math]t[/math] se encuentra en el denominador y cuando la [math]x=0[/math] la exponencial se convierte en 1 (lo cual hace que no sea capaz de corregirse la singularidad gracias a la exponencial, tal y como ocurre para otros casos de x en los que la exponencial tiende a 0 cuando t tiende a 0), provocando esa singularidad.

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);

%Dibuja las gráficas
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

15 .Referencias

• Series y transformadas de Fourier. Domingo Alcaraz Candela.
• [1]
• Serie compleja de Fourier
• [2]