Revisión del 20:39 6 mar 2024

1 Introducción

El primer objetivo de este trabajo es estudiar la ecuación del calor en un dominio acotado. Veremos también cómo afecta la variación de cada uno de los parámetros en la solución final del problema. Para ello, empezaremos resolviendo uno a modo de ejemplo. Posteriormente, modificaremos cada uno de los datos dejando fijo el resto y comparando las distintas soluciones obtenidas. Además, no sólo estudiaremos la solución matemática al problema sino que también trataremos de darle una interpretación física.

Por otro lado, consideraremos el caso de un dominio no acotado. Para ello, hablaremos de la solución fundamental de la ecuación del calor y veremos cómo obtener la solución de un problema de Cauchy dado a partir de esta. También veremos lo que es una solución autosemejante y hallaremos un ejemplo de este tipo de soluciones.

2 Preliminares

Antes de pasar a la resolución de diferentes ejemplos de la ecuación es necesario conocer los siguientes conceptos:

Ley de Fourier. Ley que establece que el flujo de transferencia de calor por conducción tiene como dirección la de mayor diferencia de temperatura. Además, este flujo va de las zonas más calientes a las más frías. Matemáticamente se puede modelizar como [math] q=-k\nabla T[/math], donde [math]k[/math] es lo que se conoce como conductividad térmica.

Conductividad térmica. Es una propiedad física que mide la capacidad de la conducción térmica de los materiales. En el Sistema Internacional de Unidades se mide en [math]\frac{J}{m s K}[/math]. A mayor sea esta constante, que cambia en función del material, mejor conductor del calor es.

Calor específico. Es la cantidad de calor que hay que suministrar a una unidad de masa para aumentar la temperatura una unidad de temperatura. Como se puede intuir de la definición se utilizan diferentes escalas. A mayor es el calor específico del material, más energía es necesaria para aumentar la temperatura. Se denota como [math]c[/math].

Coeficiente de difusión. Se define el coeficiente de difusión como el cociente [math]D=\frac{K}{c}[/math]. Por la definición de conductividad térmica y calor específico podemos concluir que a mayor es el coeficiente de difusión, menor es la energía para cambiar la temperatura. Como se verá más adelante esto influye directamente en la velocidad de los procesos de calentamiento o enfriamiento de un cuerpo.

Ecuación del Calor. La ecuación del calor es una ecuación diferencial en derivadas parciales que describe la distribución del calor, y por tanto las variaciones de la temperatura en una región a lo largo del transcurso del tiempo. Su expresión más general es:

donde [math] u : \Omega \subset \mathbb{R^n} \rightarrow \mathbb{R}, f : \Omega \subset \mathbb{R^n} \rightarrow \mathbb{R} [/math], [math]D[/math] es el coeficiente de difusión y [math]\Omega[/math] un abierto. Para ver su deducción ^[1].

Series de Fourier. Un desarrollo en serie de Fourier nos permite aproximar funciones mediante sumas de funciones periódicas, en nuestro caso usaremos senos y cosenos. Concretamente las utilizaremos para imponer que la solución que obtengamos de la ecuación verifique la condición inicial. Para más información ver ^[2].

3 Solución ecuación del calor

Una vez ya sabemos todo lo necesario para trabajar con la ecuación de calor, resolveremos un problema a modo de ejemplo y como comentamos en la introducción veremos cómo cambia la solución según modificamos los datos iniciales.

Consideremos una varilla metálica de longitud 1 m. Supongamos que esta se encuentra aislada por su superficie lateral, por lo que la conducción de calor sólo se produce en la dirección longitudinal. Además, sabemos que la temperatura inicial de la varilla es 0 ºC y en los extremos izquierdo y derecho se consigue mantener la temperatura a 0ºC y 1ºC respectivamente. Por último, supongamos que la conductividad térmica y calor específico son iguales a 1. Con todos estos datos podemos plantear el primer problema a resolver

3.1 Solución Estacionaria

Para resolver esta ecuación lo primero que debemos hacer es homogeneizarla. Es decir, obtener un nuevo problema mediante un cambio de variable de tal forma que las condiciones frontera sean la función nula en ambos extremos. Para ello, necesitamos hallar lo que se conoce como solución estacionaria.

Para entender qué representa esta solución, volvamos al problema anterior. Cuando dejamos pasar el tiempo y la temperatura en cada punto de la barra cambia según lo indica la ecuación del calor, después de un tiempo suficiente, la distribución de temperatura en la barra se mantiene constante. En este punto, hemos alcanzado la solución estacionaria. Escribámoslo en términos matemáticos:

Definición. Sea [math]T(x,t)[/math] solución de la ecuación del calor, llamamos solución estacionaria a [math]T^*(x):=\lim_{t \to \infty} T(x,t) [/math]. Esta verifica que [math]T_t^*(x)=0[/math]

Con estas ideas en mente podemos ya podemos plantear el siguiente problema cuya solución es la estacionaria:

Resolviéndola obtenemos que [math]T^*(x)=x[/math] cuya gráfica es:

Solución Estacionaria

% Input
Te=@(x) (x);    % Definición de la función
xx=0:10^(-3):1; % Discretización del dominio espacial

% Output
yy=Te(xx);      % Evalucación de la discretización
plot(xx,yy,'Color','b')     % Gráfico de la función
axis equal
xlim([0,1])
ylim([0,1])
colormap turbo
legend({'$T^{*}(x)=x$'},'Interpreter','latex','Location','southeast')
xlabel('$x$','Interpreter','latex')
ylabel('$T$','Interpreter','latex')

3.2 Resolución Problema Homogeneizado

Como habíamos dicho antes, nuestro objetivo era encontrar la solución estacionaria para poder hacer un cambio de variable que homogeneizara la ecuación y así poder resolver por separación de variables. Por tanto, haciendo el cambio de variable [math]u(x,t)= T(x,t)- T^*(x) [/math] obtenemos el problema

Se observa que efectivamente se trata del problema homogéneo que podemos resolver por separación de variables. Una vez hecho esto y deshaciendo el cambio obtenemos la siguiente solución:

A continuación se muestra la gráfica de la solución. En ella se aprecia claramente cómo se verifican las condiciones iniciales y de frontera. Además, se puede ver que pasado un tiempo se va aproximando a la solución estacionaria que habíamos calculado.

Solución ecuación del calor

% Input
N=33;           % Número de términos de la serie de Fourier
xx=0:10^(-2):1; % Discretización del dominio espacial
tt=0:10^(-2):1; % Discretización del dominio temporal

% Output
yy=T(xx,tt,N);  % Plot solución del problema
mesh(xx,tt,yy)
colormap turbo
xlabel('$x$','Interpreter','latex')
ylabel('$t$','Interpreter','latex')
zlabel('$T$','Interpreter','latex')

Donde la definición de la función [math]T[/math] es:

function Teval = T(x,t,N)
% Se introduce un vector fila x discretización espacial, un vector fila t
% discretización temporal y una N que da el número de términos de la serie
% de Fourier y devuelve la matriz solución de la EDP, donde el elemento aij
% hace referencia a la temperatura en el tiempo ti y espacio xj

    Teval=zeros(length(t),length(x));   % Inicialización de la matriz
    for n=1:N                           % Serie de Fourier
        Teval=Teval+exp(-n^(2)*pi^(2)*t')*sin(n.*pi.*x)*((-1)^n)*2./(n.*pi);
    end
    Teval=ones(1,length(t))'*x+Teval;   % Suma de la solución estacionaria 
    % y la del problema homogéneo para obtener la original
end

Una vez hemos resuelto la ecuación podemos abandonar el mundo puramente matemático y adentrarnos en la física. Nuestro nuevo objetivo es darle una interpretación a lo que hemos obtenido en el desarrollo anterior. Como se puede ver en la gráfica de la solución del apartado anterior, la barra comienza teniendo una temperatura uniforme e igual a 0ºC. Una vez comienza el experimento, de forma instantánea, el extremo derecho se calienta alcanzando 1ºC que se mantendrá constante a lo largo del tiempo. Es decir, esta ecuación junto con sus condiciones iniciales y frontera modelizan una barra que comienzan teniendo una temperatura uniforme y a la que se pone en contacto con una superficie en el extremo derecho que la calienta de forma instantánea, algo imposible en la realidad, y mantiene su temperatura constante a partir de ese momento. Además, en el lado izquierdo con un mecanismo similar se consigue que su temperatura sea la misma. Según avanza el experimento podemos ver cómo esta va aumentando en todos los puntos de la barra a excepción de los anteriores hasta que se llega a la solución descrita por la solución estacionaria. Esto concuerda con lo que sucede en el mundo físico, pues si [math] t\gt 0[/math] el extremo derecho de la barra está más caliente que su entorno, debido a esto habrá una transferencia de energía hacia la zona que esté más fría en forma de calor, por lo que su temperatura aumenta. Por otro lado, es lógico pensar que la distribución de temperatura pasado un tiempo sea continua, pues en la realidad no se ven cuerpos en equilibrio térmico que tengan saltos de temperatura en puntos muy próximos. Esto se debe a que la acción del paso del tiempo tiende a homogeneizarlos en términos de temperatura. Esta idea concuerda con lo que se ve en la solución estacionaria; una función continua que verifica las condiciones frontera.

3.3 Interpretación del flujo en los extremos

Una vez hemos entendido analíticamente el comportamiento de la situación y que efectivamente concuerda con lo que se ve en el mundo físico, veamos que dicha interpretación también podemos encontrarla de forma implícita en la solución de la ecuación. Según la ley de Fourier, dada una función temperatura de un cuerpo [math]T(x,t)[/math] suponiendo que el flujo de calor que lo atraviesa sólo puede tener una dirección, se puede calcular como:

Así, podemos hallar el flujo a lo largo de la barra:

Evaluando en los extremos obtenemos:

Para poder estudiar cómo se comporta el flujo es más fácil verlo a partir de las gráficas de las anteriores funciones que se muestran a continuación.

Flujo en los extremos

clear all
% Input
N=33;               % Número de términos de la serie de Fourier
tt=0:10^(-2):1;     % Discretización del dominio temporal

% Output
subplot(1,2,1)      % Plot gráfica flujo izquierda
plot(tt,(fi(tt,N)),'Color','blue')
xlim([0,1])
ylim([-2,0])
xlabel('$t$','Interpreter','latex')
ylabel('$q$','Interpreter','latex')
subtitle('Extremo izquierdo','Interpreter','latex')

subplot(1,2,2)      % Plot gráfica flujo derecha
plot(tt,(fd(tt,N)),'Color','blue')
xlim([0,1])
ylim([-2,0])
xlabel('$t$','Interpreter','latex')
ylabel('$q$','Interpreter','latex')
subtitle('Extremo derecho','Interpreter','latex')

Donde la función [math]fi[/math] viene dada por:

function Teval = fi(t,N)
% Se introduce un vector fila t discretización temporal y una N que da el 
% número de términos de la serie de Fourier y devuelve el vector flujo en 
% el extremo izquierdo, donde el elemento i hace referencia al instante ti

    Teval=zeros(1,length(t));           % Inicialización del vector
    for n=1:N                           % Serie de Fourier
        Teval=Teval+exp(-n^(2)*pi^(2)*t')*cos(n.*pi.*0)*((-1)^n)*2;
    end                         
    Teval=-(ones(1,length(t))'+Teval);  % Suma de la derivada de la 
    % solución estacionaria y la del problema homogéneo para obtener la 
    % original, cambiado de signo porque es el flujo
end

Y la función [math]fd[/math] viene dada por:

function Teval = fd(t,N)
% Se introduce un vector fila t discretización temporal y una N que da el 
% número de términos de la serie de Fourier y devuelve el vector flujo en 
% el extremo derecho, donde el elemento i hace referencia al instante ti

    Teval=zeros(1,length(t));           % Inicialización del vector
    for n=1:N                           % Serie de Fourier
        Teval=Teval+exp(-n^(2)*pi^(2)*t')*cos(n.*pi.*1)*((-1)^n)*2;
    end
    Teval=-(ones(1,length(t))'+Teval);  % Suma de la derivada de la 
    % solución estacionaria y la del problema homogéneo para obtener la 
    % original, cambiado de signo porque es el flujo
end

En ambas gráficas se puede ver que el flujo es siempre negativo, por lo que el calor va en el sentido negativo del eje [math]X[/math], es decir, el calor va del extremo derecho de la barra al izquierdo. Por tanto, esto coincide con el sentido físico que comentábamos en el apartado anterior. Otro hecho relevante es que la gráfica correspondiente al flujo en el extremo derecho para [math]t=0[/math] hay una asíntota vertical. Esto se debe a la discontinuidad en la función temperatura en esta frontera (recordemos [math]T(x,0)=0[/math] y [math]T(1,t)=1[/math] si [math]t\gt0[/math]) ya que para que la temperatura cambie de forma instantánea necesitaríamos "una cantidad infinita de energía". Posteriormente, el flujo se va reduciendo con el tiempo. La interpretación física de esta moderación es que al principio encontramos una gran diferencia de temperatura en un entorno del extremo derecho. Por tanto, por la Ley de Fourier debe de haber un flujo de calor de derecha a izquierda. Esto provoca un aumento en la temperatura en dichos puntos y con ello la diferencia mencionada va disminuyendo. Como esta es cada vez menor, el flujo entrante es menor también. En cuanto al extremo izquierdo podemos ver que en el instante inicial el flujo es nulo, ya que en un entorno de este punto la temperatura es constante. Según avanza el tiempo, debido a esa entrada de calor por el extremo derecho la temperatura en dicho entorno deja de ser constante, lo que provoca el flujo de calor de la zona más caliente a la más fría. Como dicha diferencia es cada vez mayor, el flujo aumenta con el tiempo hasta que se estabiliza en [math]q(0,t)=-1[/math]. Es importante recalcar que este es precisamente el mismo valor al que tiende el flujo en el otro extremo de la barra. Esto se debe a que se llega a la solución estacionaria previamente calculada, en la que la temperatura no es constante a lo largo del intervalo. Por un lado, al ser la solución estacionaria no puede haber variación de temperatura, por lo que el flujo debe ser igual en ambos extremos. Si hubiera una diferencia neta entre ambos valores, habría al menos un punto en la barra que es una fuente o sumidero de calor; se enfría o calienta respectivamente. Por otro lado, al ser una solución estacionaria en el que la temperatura no es constante a lo largo de la barra debe de existir un flujo no nulo en los extremos. Si este fuera nulo, tendríamos una barra aislada y su temperatura tendería a igualarse.

4 Efecto modificación coeficiente de difusión

Lo que hemos visto hasta ahora es la resolución matemática del problema Cauchy de un ejemplo en concreto y la comprensión física del fenómeno que modeliza. Sin embargo, como ya adelantamos en la primera sección del documento, otro de nuestros objetivos era estudiar el efecto que tenían la variación de los distintos parámetros en la solución final del problema. En este caso veremos cómo cambia esta según lo hace el coeficiente de difusión.

Como vimos en preliminares el coeficiente de difusión depende del calor específico de forma inversamente proporcional y la conductividad térmica de forma directamente proporcional. En el ejemplo anterior teníamos [math]K=1[/math] y [math]c=1[/math], por lo que [math]D=1[/math]. Supongamos ahora que la conductividad térmica se ve modificada a [math]K=\frac{1}{2}[/math] y por tanto [math]D=\frac{1}{2}[/math]. Sin embargo, en vez de trabajar directamente con este valor, hagámoslo de forma general y luego particularicemos. Así el problema general queda:

Hagamos el cambio de variable [math]\tau=\alpha t [/math] con [math]\alpha \in \mathbb{R} [/math] para simplificar lo anterior. De esta manera buscamos que [math]T(t,x)=u(\tau,x)=u(\alpha t,x)[/math], y por tanto obtenemos que :

Finalmente, tomando [math]\alpha=D[/math] llegamos a [math]u_t-u_{xx}=0[/math], que resulta la misma ecuación que ya hemos resuelto en apartados anteriores. Además, es fácil comprobar que las condiciones iniciales y frontera no cambian respecto al problema original, pues sólo estamos cambiando la escala temporal. Así, llegamos al mismo problema que teníamos antes y cuya solución conocemos:

Ahora sólo queda deshacer el cambio de variable para [math]\alpha=D[/math] obtenemos que:

Finalmente particularizamos para [math]D=\frac{1}{2}[/math]:

Para comparar con la solución anterior es mejor dibujarla, obteniendo lo que se obtiene a continuación:

Solución con coeficiente de difusión [math]D=\frac{1}{2}[/math]

% Input
% Input
N=33;           % Número de términos de la serie de Fourier
xx=0:10^(-2):1; % Discretización del dominio espacial
tt=0:10^(-2):1; % Discretización del dominio temporal

% Output
yy=Tp(xx,tt,N);  % Plot solución del problema
mesh(xx,tt,yy)
colormap turbo
xlabel('$x$','Interpreter','latex')
ylabel('$t$','Interpreter','latex')
zlabel('$T$','Interpreter','latex')

Donde la definición de la función [math]Tp[/math] es la siguiente:

function Teval = Tp(x,t,N)
% Se introduce un vector fila x discretización espacial, un vector fila t
% discretización temporal y una N que da el número de términos de la serie
% de Fourier y devuelve la matriz solución de la EDP, donde el elemento aij
% hace referencia a la temperatura en el tiempo ti y espacio xj

    Teval=zeros(length(t),length(x));   % Inicialización de la matriz
    for n=1:N                           % Serie de Fourier
        Teval=Teval+exp(-n^(2)*pi^(2)*t'/2)*sin(n.*pi.*x)*((-1)^n)*2./(n.*pi);
    end
    Teval=ones(1,length(t))'*x+Teval;   % Suma de la solución estacionaria 
    % y la del problema homogéneo para obtener la original
end

En esta no se ve un cambio aparente respecto de la anterior gráfica. Sin embargo, podemos apreciar mejor este cambio si para un punto fijado vemos la evolución de la diferencia entre el estado estacionario y las soluciones obtenidas para cada uno de los valores de la difusión. De esta manera podemos comparar la velocidad de calentamiento de la barra en un punto; hagámoslo para [math]x=\frac{1}{2}[/math].

Evolución de las diferencias de temperatura entre el estado estacionario y las soluciones para [math]D=1[/math] y [math]D=\frac{1}{2}[/math]

% Input
N=33;                                   % Número de términos de la serie de Fourier
tt=0:10^(-2):1;                         % Discretización del dominio temporal

% Output
subplot(1,2,1)                          % Plot problema original
yy=1/2*ones(length(tt),1)-T(1/2,tt,N);  % Evaluación de la diferencia
plot(tt,yy,'Color','blue')                             % Plot de la función
xlabel('$t$','Interpreter','latex')
ylabel('$T$','Interpreter','latex')
subtitle('Problema original','Interpreter','latex')

subplot(1,2,2)                          % Plot problema nuevo
yy=1/2*ones(length(tt),1)-Tp(1/2,tt,N);  % Evaluación de la diferencia
plot(tt,yy,'Color','blue')                             % Plot de la función
xlabel('$t$','Interpreter','latex')
ylabel('$T$','Interpreter','latex')
subtitle('Problema modificado','Interpreter','latex')

Ahora sí podemos ver que para un coeficiente de difusión mayor la diferencia de las temperaturas se acerca más rápido a 0 que en la otra ,es decir, que tiende más rápido a la solución estacionaria. Por tanto, se puede afirmar que la diferencia entre ambas situaciones radica en la velocidad del proceso. Esto cuadra con nuestra intuición física del problema, pues como se ha mencionado en preliminares a mayor conductividad térmica más fácil es para el calor atravesar el material del que está hecho la barra y por tanto más rápido es el proceso descrito. Como dicho coeficiente ha disminuido a la mitad, el material con el que está hecha la segunda barra es peor conductor del calor y por tanto el proceso es más lento que en el primer caso, como se ha podido observar en la gráfica.

5 Efecto modificación de la condición inicial

Tras visualizar el efecto que tiene la modificación del coeficiente de difusión, veamos el efecto de un cambio en la condición inicial. Tomemos ahora el siguiente problema.

La condición inicial por tanto ahora es esta función:

Gráfica de la nueva condición inicial

u0=@(x) (max(0,1-4*abs(x-1/2)));   % Definición de la función
xx=0:10^(-3):1;                 % Discretización del dominio espacial
yy=u0(xx);                      % Evalucación de la discretización
plot(xx,yy,'Color','blue')                     % Gráfico de la función
axis equal
colormap turbo
legend({'$u_{0}(x)$'},'Interpreter','latex','Location','southeast')
xlabel('$x$','Interpreter','latex')
ylabel('$T$','Interpreter','latex')

Esta representa un estado inicial de la barra con mayor temperatura entorno a su centro. En este caso las condiciones frontera son nulas por lo que no hace falta homogeneizar el problema como hemos hecho anteriormente. Buscamos entonces una solución por el método de separación de variables y llegamos a

Aquí [math]c_{k}[/math] denotan los coeficientes de Fourier de la función [math]u(x,0) = \max\lbrace 0,1−4\lvert x−\frac{1}{2}\rvert\rbrace[/math] en la base [math]\lbrace sen(k\pi x)\rbrace_{n\in\mathbb{N}}[/math] que calcularemos mediante la siguiente integral aproximándola numéricamente en MatLab:

Para más información ver ^[3]. La gráfica de la solución quedaría de la siguiente manera:

Gráfica de la solución

% Input
N=33;           % Número de términos de la serie de Fourier
xx=0:10^(-2):1; % Discretización del dominio espacial
tt=0:10^(-2):1; % Discretización del dominio temporal

% Output
yy=T2(xx,tt,N); % Plot solución del problema
mesh(xx,tt,yy)
colormap turbo
xlabel('$x$','Interpreter','latex')
ylabel('$t$','Interpreter','latex')
zlabel('$T$','Interpreter','latex')

A continuación se muestra la definición de la función [math]T2[/math]

function Teval = T2(x,t,N)
% Se introduce un vector fila x discretización espacial, un vector fila t
% discretización temporal y una N que da el número de términos de la serie
% de Fourier y devuelve la matriz solución de la EDP, donde el elemento aij
% hace referencia a la temperatura en el tiempo ti y espacio xj

    Teval=zeros(length(t),length(x));   % Inicialización de la matriz
    for n=1:N                           % Serie de Fourier
        Teval=Teval+trapz(max(0,1-4*abs(x-1/2)).*sin(pi*n*x))*exp(-n^(2)*pi^(2)*t')*sin(n.*pi.*x)/trapz(sin(pi*n*x).^(2));
    end
end

Podemos ver en esta gráfica cómo en los primeros instantes el calor está concentrado en el centro de la barra y cómo se distribuye a lo largo de esta. Por tanto, el máximo disminuye con el tiempo. Debido a que la temperatura en el interior de la barra es positiva y en sus extremos nula, por la Ley de Fourier llega a un estado estacionario de temperatura constante igual a cero. Veamos el flujo en los extremos para comprender mejor cómo se disipa el calor.

Flujo en los extremos en la barra en función del tiempo

clear all
% Input
N=33;                   % Número de términos de la serie de Fourier
xx=0:10^(-2):1;         % Discretización del dominio espacial
tt=0:10^(-2):1;         % Discretización del dominio temporal

%Output
subplot(1,2,1)          % Plot gráfica flujo izquierda
plot(tt,(fi2(xx,tt,N)),'Color','blue')
xlim([0,1])
ylim([-0.2,0.2])
xlabel('$t$','Interpreter','latex')
ylabel('$q$','Interpreter','latex')
subtitle('Extremo izquierdo','Interpreter','latex')

subplot(1,2,2)          % Plot gráfica flujo derecha
plot(tt,(fd2(xx,tt,N)),'Color','blue')
xlim([0,1])
ylim([-0.2,0.2])
xlabel('$t$','Interpreter','latex')
ylabel('$q$','Interpreter','latex')
subtitle('Extremo derecho','Interpreter','latex')

Se puede ver cómo en la izquierda el flujo comienza a ser negativo, lo que implica que el calor se desplaza hacia la izquierda en dicho punto, y en la derecha comienza siendo positivo, lo que implica que en dicho punto se mueve hacia la derecha. Esto nos describe de nuevo el fenómeno enunciado de cómo el calor se va escapando por los extremos. A medida que se recorre el eje temporal en ambas gráficas se ve como cada vez se escapa menos calor hasta que el flujo es cero; ni entra ni sale calor. Este comportamiento hace referencia a cuando se alcanza el estado estacionario de temperatura igual a cero y no hay transferencia de energía en ningún sentido.

6 Efecto modificación condiciones frontera

Tras la modificación de la condición inicial, veamos como influyen las condiciones frontera de distinto tipo. Pongamos ahora una condición del tipo [math] T_x=0[/math] en el extremo derecho. La interpretación física de esta condición sería la de aislar el extremo derecho para que no entre ni salga calor a través de este, ya que el flujo es nulo.

Resolviendo de nuevo por separación de variables obtenemos que la solución es la función

Los [math]c_{k}[/math] denotan de nuevo los coeficientes de Fourier de la función [math]T(x,0) = \max\lbrace 0,1−4\lvert x−\frac{1}{2}\rvert\rbrace[/math] en la base [math]\lbrace sen(\frac{\pi}{2}(2n-1) x)\rbrace_{n\in\mathbb{N}}[/math]. Grafiquemos la solución para ver su desarrollo en tiempo y espacio:

Gráfica de la solución

% Input
N=33;           % Número de términos de la serie de Fourier
xx=0:10^(-2):1; % Discretización del dominio espacial
tt=0:10^(-2):1; % Discretización del dominio temporal

% Output
yy=T3(xx,tt,N); % Plot solución del problema
mesh(xx,tt,yy)
colormap turbo
xlabel('$x$','Interpreter','latex')
ylabel('$t$','Interpreter','latex')
zlabel('$T$','Interpreter','latex')

A continuación se muestra la definición de [math]T3[/math]:

function Teval = T3(x,t,N)
% Se introduce un vector fila x discretización espacial, un vector fila t
% discretización temporal y una N que da el número de términos de la serie
% de Fourier y devuelve la matriz solución de la EDP, donde el elemento ij
% hace referencia a la temperatura en el tiempo ti y espacio xj

    Teval=zeros(length(t),length(x));   % Inicialización de la matriz
    for n=1:N                           % Serie de Fourier
        Teval=Teval+trapz(max(0,1-4*abs(x-1/2)).*sin(pi*(2*n-1)*x/2))*exp(-pi^(2)*(2*n-1)^(2)*t'/4)*sin(pi*(2*n-1)*x/2)/trapz(sin(pi*(2*n-1)*x/2).^(2));
    end
end

Llama la atención cómo en los instantes iniciales es igual que en el problema anterior. El máximo va disminuyendo mientras la temperatura en el resto de la barra aumenta. En el extremo izquierdo se observa una tendencia similar al problema previo, mientras que en el extremo derecho aumenta la temperatura ya que está aislado. Pasado un tiempo la barra acaba tendiendo de igual manera al estado estacionario de temperatura cero, con la diferencia de que el calor se va perdiendo solo por uno de los dos extremos. Veamos un gráfico del flujo en los extremos para observar mejor esta tendencia.

Flujo en los extremos en función del tiempo

clear all
% Input
N=33;                   % Número de términos de la serie de Fourier
xx=0:10^(-2):1;         % Discretización del dominio espacial
tt=0:10^(-2):1.5;         % Discretización del dominio temporal

%Output
subplot(1,2,1)          % Plot gráfica flujo izquierda
plot(tt,(fi3(xx,tt,N)),'Color','blue')
xlim([0,1.5])
ylim([-1,0.2])
xlabel('$t$','Interpreter','latex')
ylabel('$q$','Interpreter','latex')
subtitle('Extremo izquierdo','Interpreter','latex')

subplot(1,2,2)          % Plot gráfica flujo derecha
plot(tt,(fd3(xx,tt,N)),'Color','blue')
xlim([0,1.5])
ylim([-1,0.2])
xlabel('$t$','Interpreter','latex')
ylabel('$q$','Interpreter','latex')
subtitle('Extremo derecho','Interpreter','latex')

A continuación se muestran las definiciones de las funciones [math]fi3[/math] y [math]fd3[/math].

function Teval = fi3(x,t,N)
% Se introduce un vector fila t discretización temporal y una N que da el 
% número de términos de la serie de Fourier y devuelve el vector flujo en 
% el extremo izquierdo, donde el elemento i hace referencia al instante ti

    Teval=zeros(length(t),1);       % Inicialización del vector
    for n=1:N                       % Serie de Fourier
        Teval=Teval+pi*(2*n-1)*trapz(max(0,1-4*abs(x-1/2)).*sin(pi*(2*n-1)*x/2))*exp(-pi^(2)*(2*n-1)^(2)*t'/4)*cos(pi*(2*n-1)*0/2)/(2*trapz(sin(pi*(2*n-1)*x/2).^(2)));
    end
    Teval=-Teval;                   % Cambio de signo porque es un flujo
end

function Teval = fd3(x,t,N)
% Se introduce un vector fila t discretización temporal y una N que da el 
% número de términos de la serie de Fourier y devuelve el vector flujo en 
% el extremo derecho, donde el elemento i hace referencia al instante ti

    Teval=zeros(length(t),1);       % Inicialización del vector
    for n=1:N                       % Serie de Fourier
        Teval=Teval+pi*(2*n-1)*trapz(max(0,1-4*abs(x-1/2)).*sin(pi*(2*n-1)*x/2))*exp(-pi^(2)*(2*n-1)^(2)*t'/4)*cos(pi*(2*n-1)*1/2)/(2*trapz(sin(pi*(2*n-1)*x/2).^(2)));
    end
    Teval=-Teval;                   % Cambio de signo porque es un flujo
end

Se puede observar en primer lugar el flujo nulo del extremo derecho impuesto por la condición frontera; como no podía ser de otra manera. En el extremo izquierdo vemos cómo se va escapando todo el calor; el flujo es negativo. Comparándolo con el ejemplo anterior tarda más en anularse ya que todo el calor de la barra sale por ese extremo y no por ambos a la vez. Finalmente, se puede observar cómo este acaba haciéndose nulo, y dando por resultado que toda la barra se encuentre a cero grados, el estado estacionario.

7 Principio del máximo

Como se ha visto en los gráficos anteriores, la ecuación del calor tiene un efecto regularizante que hace que la función en cuanto se despega del instante inicial pasa a ser infinitamente derivable. Esto intuitivamente nos indica que los picos de temperatura se deberían encontrar en el instante inicial de la barra o a lo largo de los extremos de esta, la llamada frontera parabólica. Formalmente, esta intuición es un resultado conocido como el principio del máximo.

Teorema (Principio del máximo). Sea [math]T\in\mathcal{C}^{2,1}(\mathcal{Q}_T)\times\mathcal{C}(\overline{\mathcal{Q}_T})[/math] que cumple

entonces

donde [math]\partial_{p}\mathcal{Q}_T[/math] denota la frontera hiperbólica, conformada por la frontera de [math]{Q}_T[/math] excluyendo los puntos de la forma [math](x,T)[/math].

Rescatemos las soluciones de los apartados anteriores para comprobar que se verifica este teorema.

PRIMERA SOLUCION En la primera solución se puede observar que la imposición de forzar la temperatura del extremo derecho a que sea constantemente igual a [math]1[/math] hace que el máximo se encuentre en este. El mínimo lo podemos ver en el otro extremo, ya que estamos forzando a que se mantenga a temperatura cero. Adicionalmente se ve que esta temperatura también se alcanza en el instante inicial. Tanto el máximo como el mínimo se alcanzan en la frontera parabólica, como nos dicta el principio del máximo.

SEGUNDA SOLUCION En la gráfica de la segunda solución se puede observar que el máximo se alcanza esta vez en el instante inicial, ya que a medida que avanza el tiempo el calor se va disipando hacia los extremos. La imposición de que los extremos estén constantemente a temperatura igual a cero tiene como consecuencia que el mínimo se halle sobre estos. De manera añadida también se puede encontrar este valor en la condición inicial. De nuevo, los casos extremos se observan en la frontera parabólica.

TERCERA SOLUCION La tercera solución es bastante similar a la anterior, por ello el valor máximo se halla de igual manera en el instante inicial pese a un amago de máximo en una acumulación de calor en el extremo derecho, fruto de la condición frontera de flujo nulo. De nuevo, el mínimo se encuentra en la condición inicial y en el extremo derecho. Transcurrido un tiempo lo podemos hallar también en el extremo derecho de la barra. De nuevo todos estos casos se observan en la frontera parabólica.

8 Solución Fundamental de la ecuación del calor

En los apartados anteriores siempre que hemos resuelto la ecuación del calor lo hemos hecho imponiendo unas condiciones iniciales y de frontera, por lo que el dominio de definición de la solución está acotado. ¿Se podría encontrar una solución sin imponer alguna de esas condiciones? La respuesta a esta pregunta es afirmativa, se puede resolver sin imponer condiciones de frontera; la solución no está acotada en la variable espacial. A esta función se le conoce como solución fundamental de la ecuación del calor por una propiedad que tiene que estudiaremos más adelante. La expresión de dicha solución, para dimensión [math]n=1[/math] y coeficiente de difusión [math]D=1[/math] es la que sigue:

Para ver cómo se obtiene ^[4]. A continuación se muestra la gráfica de esta función.

Gráfica de la solución fundamental de la ecuación del calor

% Input
xx=-1:10^(-1.8):1;                                  % Discretización del espacio
tt=10^(-2):10^(-2.7):1;                             % Discretización del tiempo

% Output
u=@(x,t)(exp(-x.^(2)./(4*t'))./(2*sqrt(pi*t')));    % Definición de la solución fundamental
yy=u(xx,tt);                                        % Evaluación de la solución fundamental
mesh(xx,tt,yy)                                      % Grafica de la solución
colormap turbo
xlabel('$x$','Interpreter','latex')
ylabel('$t$','Interpreter','latex')

9 Soluciones Autosemejantes

Otra tipo de soluciones de especial interés son las autosemejantes, estas son las que son invariantes bajo ciertos cambios de escala en espacio y tiempo. Un tipo de estas soluciones autosemejantes son de la forma [math]u(x,t)=U(\frac{x}{\sqrt{t}})[/math]. Veamos ahora una expresión más precisa que describe este tipo de funciones:

Como [math] U(\frac{x}{\sqrt{t}}) [/math] es solución de la ecuación, derivando llegamos a la siguiente expresión:

Haciendo el cambio de variable [math]\xi=\frac{x}{\sqrt{t}}[/math] llegamos a:

Operando un poco llegamos a la EDO:

Resolviéndola y deshaciendo el cambio de variable llegamos a la expresión buscada:

Nótese que al igual que para la solución fundamental no hemos necesitado imponer unas condiciones frontera. Es más, no hemos empleado tampoco ninguna condición inicial. Intentemos dibujar esta función para ver su comportamiento. Para ello, debemos imponer unas condiciones que nos permitan hallar los valores de las constantes, impongamos lo siguiente y veamos qué obtenemos:

Es fácil ver que la segunda y tercera expresión resultan en que [math]c=0[/math]. Por otro lado, a partir de la primera tenemos que [math]A+2c=1[/math] y por tanto [math]A=1[/math]. Sustituyendo estos valores en la solución que hemos obtenido anteriormente obtenemos la expresión que buscábamos:

Finalmente podemos dibujar su gráfica:

Solución autosemejante de la ecuación del calor

% Input
xx=0:10^(-2):5;                     % Discretización del espacio
tt=10^(-3):10^(-2.7):1;             % Discretización del tiempo

% Output
u=@(t,x) exp(-x./((4*t').^(1/2)));  % Definición de la solución
yy=u(tt,xx);                        % Evaluación de la solución
mesh(xx,tt,yy)                      % Grafica de la solución
colormap turbo
xlabel('x','Interpreter','latex')
ylabel('t','Interpreter','latex')

10 Cálculo de soluciones a partir de la solución fundamental

Término de la base [math] \cos(n\pi x)[/math]

% Input
xx=-10:10^(-1):10;                                  % Discretización del espacio
tt=10^(-3):10^(-2.7):10;                            % Discretización del tiempo
% Output
K=@(x,t)(exp(-x.^(2)./(4*t'))./(2*sqrt(pi*t')));    % Definición de la solución fundamental
ff=u3(xx,tt);                                       % Cálculo de la solución
mesh(xx,tt,ff)                                      % Gráfica de la solución
colormap turbo
xlabel('$x$','Interpreter','latex')
ylabel('$t$','Interpreter','latex')

function Teval = u3(x,t)
% Se introduce un vector x discretización especial y un vector t
% discretización temporal y se devuelve la solución del problema con
% condición inicial u0 en forma de matriz, donde el elemento ij hace
% referencia a la temperatura en el espacio xj en el instante de tiempo ti
    paso=10^(-1);                                       % Separación de la discretización para la integral
    y=-1:paso:1;                                        % Discretización para la integral
    ly=length(y);
    K=@(x,t)(exp(-x.^(2)./(4*t'))./(2*sqrt(pi*t')));    % Definición de la solución fundamental
    Teval=zeros(length(t),length(x));                   % Inicialización de la matriz
    for i=1:length(x)                                   % Evaluación de la integral para cada x y cada t
        for j=1:length(t)
            Teval(j,i)=trapz(K(x(i)*ones(1,ly)-y,t(j)).*u0(y));
        end
    end
end

function Teval = u0(x)
% Se introduce un vector fila x y se devulve su evaluación discreta
% por la función u0
    Teval=zeros(1,length(x));   % Inicialización del vector evaluado
    for i=1:length(x)           % Evaluación inicial de cada componente
        if abs(x(i))<=1
            Teval(i)=1;
        end
    end
end

Término de la base [math] \cos(n\pi x)[/math]

% Input
xx=-5:10^(-1):5;                                    % Discretización del espacio
% Output
K=@(x,t)(exp(-x.^(2)./(4*t'))./(2*sqrt(pi*t')));    % Definición de la solución fundamental
ff=u3(xx,tt);

subplot(1,3,1)                                      % t=10^(-3)
plot(xx,u3(xx,10^(-3)),'Color','blue')              % Gráfica de la solución
xlabel('$x$','Interpreter','latex')
xlim([xx(1),xx(end)])
ylim([0,12])
subtitle('$t=10^{-3}$','Interpreter','latex')

subplot(1,3,2)                                      % t=10^(-2)
plot(xx,u3(xx,10^(-2)),'Color','blue')              % Gráfica de la solución
xlabel('$x$','Interpreter','latex')
xlim([xx(1),xx(end)])
ylim([0,12])
subtitle('$t=10^{-2}$','Interpreter','latex')

subplot(1,3,3)                                      % t=10^(-1)
plot(xx,u3(xx,10^(-1)),'Color','blue')              % Gráfica de la solución
xlabel('$x$','Interpreter','latex')
xlim([xx(1),xx(end)])
ylim([0,12])
subtitle('$t=10^{-1}$','Interpreter','latex')

11 Solución fundamental dimensión 2

Término de la base [math] \cos(n\pi x)[/math]

% Input
xx=-1:10^(-1.8):1;                                  % Discretización del espacio
yy=-1:10^(-1.8):1;                                  % Discretización del tiempo

% Output
subplot(1,3,1)                                      % t=10^(-3)
zz=K2(xx,yy,10^(-3));                               % Evaluación de la solución fundamental
mesh(xx,yy,zz)                                      % Gráfica de la solución
colormap turbo
xlabel('$x$','Interpreter','latex')
ylabel('$y$','Interpreter','latex')
subtitle('$t=10^{-3}$','Interpreter','latex')

subplot(1,3,2)                                      % t=10^(-2)
zz=K2(xx,yy,10^(-2));                               % Evaluación de la solución fundamental
mesh(xx,yy,zz)                                      % Gráfica de la solución
colormap turbo
xlabel('$x$','Interpreter','latex')
ylabel('$y$','Interpreter','latex')
subtitle('$t=10^{-2}$','Interpreter','latex')

subplot(1,3,3)                                      % t=10^(-1)
zz=K2(xx,yy,10^(-1));                               % Evaluación de la solución fundamental
mesh(xx,yy,zz)                                      % Gráfica de la solución
colormap turbo
xlabel('$x$','Interpreter','latex')
ylabel('$y$','Interpreter','latex')
subtitle('$t=10^{-1}$','Interpreter','latex')

Término de la base [math] \cos(n\pi x)[/math]

function Teval = K2(x,y,t)
% Se introduce un vector fila x discretización espacial, un vector fila y
% un instante de tiempo t y se devuelve la solución fundamental bidimensional
% en el instante de matriz, donde el elemento ij hace referencia a la temperatura 
% en el espacio xj,yi
    Teval=zeros(length(y),length(x));                   % Inicialización de la matriz
    for i=1:length(x)                                   % Evaluación de la solución para cada x, y, t
        for j=1:length(y)
            Teval(j,i)=(exp(-(x(i).^(2)+y(j).^(2))./(4*t))./(4*pi*t));
        end
    end
end

11.1 Referencias

1. ↑ Salsa, S., Verzini, G. (2022). Partial Differential Equations in Action: From Modelling to Theory. Alemania: Springer International Publishing.
2. ↑ [1] Gomez, C., Rodriguez, J. (2024). Series de Fourier
3. ↑ [2] Gomez, C., Rodriguez, J. (2024). Series de Fourier
4. ↑ Salsa, S., Verzini, G. (2022). Partial Differential Equations in Action: From Modelling to Theory. Alemania: Springer International Publishing. Pg 36-38

[.[Categoría:Grado en Matemáticas].] [.[Categoría:EDP].] [.[Categoría:EDP23/24].]