Revisión del 18:06 6 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á:

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.

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] que verifica:

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

Nuestra solución correspondiente al apartado anterior cumplen 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] además por el apartado anterior sabemos que la u que es la solución correspondiente al sistema cumple la expresión [math]u_t - \Delta u\leq 0 [/math]. Además la condición inicial converge uniformemente de manera que la condición que necesitamos para la continuidad se cumple a la perfección. (HAY QUE JUSTIFICAR POR QUE LA CONDICION INICIAL CONVERGE UNIFORMEMENTE).

De esta forma si consideramos dos soluciones u y v de tal forma que definimos w=u-v que satisface:

Aplicando el principio del máximo 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 aqui deducimos que [math]w \equiv 0[/math] por tanto podemos concluir que u=v, 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, para poder hacer una buena representación tomaremos [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:

11.1 .Código

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, dichas condiciones serán; la condición inicial [math] u(x, 0) = 0 [/math] y la condición frontera [math]u(0, t) = 1 [/math]. La solución de [math]u_{t} −u_{xx} = 0[/math] en [math] x\gt0, t\gt0[/math] satisfaciendo las condiciones u(0,t)=1, cuando [math]u(x,t) \rightarrow 0 [/math] si [math]x \rightarrow \infty[/math] para [math] x\gt0 [/math] y [math] u(x,0)=0[/math] para [math] x\gt0 [/math].

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á [math]u(x,t)=U(\frac{x}{\sqrt{t}})[/math]. Sustituyendo esa expresión en la ecuación del calor obtenemos lo siguiente:

​ Se trata de una EDO lineal de orden 2 cuya solución será de la forma:

Donde si hallamos tanto [math]\lambda_{1}[/math] como [math]\lambda_{2}[/math] obtenemos los siguientes valores:

Luego la solución final que obtenemos será:

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.

13.1 .Código

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