Diferencia entre revisiones de «Ecuación del calor, equipo LUA»

Revisión del 22:17 7 mar 2024

Trabajo realizado por estudiantes
Título	Ecuación del calor. Grupo LUA
Asignatura	EDP
Curso	2023-24
Autores	Luis Carreras Hoyos, Lucía Gil Ruiz y Alejandra Hernández Sieber
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

Contenido

1 Introducción
2 Conceptos previos
3 Modelización de la ecuación del calor con una dimensión
4 Cambio de coeficientes
- 4.1 Simulación
5 Cambio de condiciones
6 Cambios de condiciones de frontera
- 6.1 Simulación
- 6.2 Principio del máximo
7 Solución fundamental de la ecuación del calor
8 Referencias

1 Introducción

La ecuación del calor es una de las ecuaciones diferenciales parciales más fundamentales en la física y la ingeniería, utilizada para modelar la difusión del calor en sólidos, líquidos y gases. Su estudio proporciona una comprensión profunda de los procesos de transferencia de calor y la evolución temporal de la temperatura en sistemas físicos. En este trabajo, nos enfocamos inicialmente en la ecuación del calor en una dimensión, explorando su formulación matemática y calculando soluciones para diversas condiciones iniciales y de frontera, además de examinar las implicaciones físicas de estas soluciones. En concreto, estudiaremos la variación temporal de la temperatura en función de la posición a lo largo de una varilla metálica, considerando condiciones iniciales y de frontera diversas, como condiciones de contorno de tipo Dirichlet y Neumann. Posteriormente, ampliaremos nuestro enfoque para abordar la ecuación del calor en dos dimensiones. Aunque este estudio será menos extenso en comparación con el caso unidimensional, nos permitirá explorar cómo se extienden los conceptos y métodos desarrollados previamente a sistemas bidimensionales.

2 Conceptos previos

Antes de comenzar con el trabajo, en esta sección vamos a definir algunos conceptos previos necesarios para la completa comprensión del mismo.

2.1 Contextualización y deducción de la ecuación del calor

Procedemos a deducir la ecuación del calor. En primer lugar, la Ley de enfriamiento de Newton lleva el nombre de Isaac Newton, quien la formuló en el siglo XVII. Esta es una ley empírica cuya idea es que la velocidad a la que disminuye (o aumenta) la temperatura [math] T(t) [/math] de un cuerpo es proporcional a la diferencia entre la temperatura ambiental [math] T_a [/math] y la del cuerpo. Es decir,

[math] \frac{dT(t)}{dt} = K \cdot (T(t)-T_a) ,~~ T(0)=T_0, ~~~ K \lt0 [/math]

Posteriormente, J.Fourier estableció que dichas temperaturas no eran constantes. Consideró [math] V [/math] un volumen de control y definió [math] e(\mathbf{x},t) [/math] como la densidad de energía calorífica, siendo esta la cantidad de energía térmica por unidad de volumen en un medio, con unidades [math] \frac{cal}{m^3} [/math].

Dado que el calor es una forma de energía, es natural hacer uso de la ley de conservación de la energía durante la deducción de la ecuación del calor. Esta postula que sobre [math] V [/math], la variación de la energía calorífica se debe al balance entre el calor que entra y el que sale, más la consideración de una posible producción externa. Por tanto, se deduce:

[math] \frac{d}{dt} \int_V e(\mathbf{x},t) \cdot d\mathbf{x} = - \int_{\partial V} \overrightarrow{q} \cdot d\overrightarrow{S} + \int_V r \cdot d\mathbf{x} [/math]

Siendo [math] \overrightarrow{q} [/math] el flujo de calor, medido en [math] \frac{cal}{m^2s} [/math] y [math] r [/math] una posible fuente externa de calor, medida en [math] \frac{cal}{m^3s} [/math].

A continuación, por el Teorema de la divergencia se tiene que sea [math] \Omega \in \mathbb{R^3} [/math] un dominio Lipschitz, entonces:

[math] \int \int \int_{\Omega} div(\overrightarrow{u} )\cdot dV = \int \int_{\partial \Omega} \overrightarrow{u} \cdot d\overrightarrow{S}[/math]

Aplicando a la Ley de conservación de la energía el Teorema de la divergencia, obtenemos

[math] \frac{d}{dt} \int_V e(\mathbf{x},t) \cdot d\mathbf{x} = - \int_V div(\overrightarrow{q}) \cdot d\mathbf{x} + \int_V r \cdot d\mathbf{x} ~~ \Longrightarrow ~~ \int_V \left[ \frac{\partial e(\mathbf{x},t)}{\partial t} + div(\overrightarrow{q}) – r \right] d\mathbf{x} = 0 ~~ \forall V \in \Omega[/math]

Como se cumple para todo volumen [math] V \in \Omega[/math], el integrando debe ser nulo, llegando a que [math] \frac{\partial e(\mathbf{x},t)}{\partial t} + div(\overrightarrow{q})=r [/math].

Para continuar, asumimos lo siguiente:

Ley de Fourier de la conducción del calor: el flujo de calor [math] \overrightarrow{q} [/math] a través de un material es directamente proporcional al producto del gradiente de temperatura [math]T(\mathbf{x},t)[/math] y la conductividad térmica [math] k [/math] del material. En términos de fórmulas, sería [math] \overrightarrow{q} = -k \cdot \nabla T(\mathbf{x},t) [/math]

La energía térmica es proporcional a la temperatura, esto es [math] e(\mathbf{x},t) = c \cdot T(\mathbf{x},t) [/math], siendo [math] c [/math] la constante de calor específico.

Con ello, hemos llegado a que

[math] \frac{\partial e(\mathbf{x},t)}{\partial t} + div(\overrightarrow{q})=r ~~ \Longrightarrow ~~ c \cdot \frac{\partial T (\mathbf{x},t)}{t} – div(k \cdot \nabla T(\mathbf{x},t)) = r ~~ \Longrightarrow ~~ c \cdot \frac{\partial T (\mathbf{x},t)}{t} – k \cdot \Delta T(\mathbf{x},t) = r[/math]

Donde en la última igualdad hemos empleado que [math] div(\nabla F) = \Delta F [/math]. Agrupando términos, se concluye que la ecuación del calor </math> sería:

[math] \frac{T(\mathbf{x},t)}{dt} - \frac{k}{c} \cdot \frac{\partial^2 T(\mathbf{x},t)}{\partial^2 x} = r [/math]

2.2 Constante de conductividad térmica y calor específico

A continuación, introducimos tres definiciones:

Constante de conductividad térmica: La constante de conductividad térmica, a la que hemos denotado como [math] k [/math], es una propiedad intrínseca de los materiales la cual describe su capacidad para conducir el calor y se mide en [math] \frac{cal}{ºC s m}[/math]. Cuanto mayor sea la constante de conductividad térmica, mayor será la capacidad del material para conducir el calor.
Calor específico: El calor específico, denotado previamente como [math] c [/math], es una propiedad de cada material que indica cuánto calor es necesario para elevar la temperatura de una unidad de masa de esa sustancia en una unidad de temperatura, medida [math] \frac{cal}{ºC g}[/math]. Cuanto mayor sea el calor específico de una sustancia, más calor se requiere para aumentar su temperatura.
Constante de difusión térmica: La constante de difusión térmica, a la que hemos denotado como [math] D [/math], es el cociente de la constante de conductividad térmica entre el calor específico.

2.3 Principio del máximo

El principio del máximo para la ecuación del calor establece que si la función [math]u(x,t) \in C^{2,1}(Q_{T}) \cap C(\overline{Q}_{T}) [/math] cumple la condición [math]u_{t}-D\Delta u \leq 0 [/math] para [math](x,t) \in Q_{T}[/math]. Entonces [math] u(x,t) [/math] verifica [math] \max_{x\in \partial_p Q_{T}} u(x,t)= \max_{x\in \overline{Q_{T}}} u(x,t) [/math].

Donde [math]Q_{T}= \Omega \times (0,T) [/math], siendo [math] \Omega \subset R^{3}[/math].

2.4 Convolución

La convolución es una operación matemática que combina dos funciones para describir la superposición entre ambas. La convolución de las funciones [math] f [/math] y [math] g [/math] se denota por [math] f *g [/math] y se define por:

[math] (f*g) (t) = \int_{-\infty}^{\infty} f(\eta) g(t-\eta) d\eta [/math]

En este trabajo, usaremos esta operación para obtener una solución de la ecuación del calor dado un dato inicial concreto, como veremos posteriormente.

3 Modelización de la ecuación del calor con una dimensión

Para estudiar el comportamiento de la ecuación del calor, vamos a comenzar tratando el siguiente problema:

Se considera una varilla metálica que ocupa el intervalo [math] [0,1] [/math] 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 [math] 0 °C [/math]. En el extremo izquierdo se consigue mantener la temperatura a [math] 0 °C[/math] mientras que en el derecho la temperatura es siempre de [math] 1°C [/math].

Aplicando la fórmula de la ecuación del calor, definiendo [math] u(x,t) [/math] como la función que determina la temperatura y tomando como constante de conductividad térmica [math] k=1~\frac{cal}{ºC s m} [/math] y calor específico [math] c=1~\frac{cal}{ºC g} [/math], modelizamos el comportamiento de la temperatura mediante el correspondiente sistema de EDP: [math] \begin{cases} u_t-u_{xx}=0 & x\in(0,1),t\gt0\\ u(x,0)=0&x\in(0,1)\\ u(0,t)=0&t\gt0\\ u(1,t)=1&t\gt0 \end{cases}[/math]

En primer lugar, cuando [math] t \rightarrow \infty[/math], [math] u_t \sim 0 [/math] y por ello, la solución no varía en tiempo, es decir, [math] u(t,x) \sim v(x) [/math], siendo [math] v(x) [/math] la solución estacionaria y cumpliendo que:

[math] \begin{cases} v_{xx}=0 & x\in(0,1)\\ v(0)=0\\ v(1)=1\\ \end{cases}[/math]

Resolviendo el sistema obtenemos [math] v(x)=x [/math], cuya representación es:

Representación de la solución estacionaria, [math] v(x)=x [/math]

clear all
close all
 
% DEFINIMOS LA SOLUCIÓN ESTACIONARIA
sol_estacionaria=@(x) x;
 
% DATOS PARA LA DISCRETIZACIÓN DEL ESPACIO
extremo_izquierdo=0;
extremo_derecho=1;
numero_puntos=100;
 
% DISCRETIZACIÓN DEL ESPACIO Y DEL TIEMPO
X=linspace(extremo_izquierdo,extremo_derecho,numero_puntos); % vector del mallado espacial
 
plot(X,sol_estacionaria(X))
xlabel('Espacio');
ylabel('Temperatura');
title("Solución estacionaria")

Esta solución estacionaria se resta a la solución original obteniendo el cambio de variable [math] w(x,t)=u(x,t)-v(x) [/math]. Mediante este cambio, obtenemos nuestro sistema de EDP’s homogeneizado:

[math] \begin{cases} w_t-w_{xx}=0 & x\in(0,1),t\gt0\\ w(x,0)=-x&x\in(0,1)\\ w(0,t)=0&t\gt0\\ w(1,t)=0&t\gt0 \end{cases}[/math]

Tomamos ahora el cambio [math] z(x,t)=-w(x,t) [/math], donde obtenemos el sistema final: [math] \begin{cases} z_t-z_{xx}=0 & x\in(0,1),t\gt0\\ z(x,0)=x&x\in(0,1)\\ z(0,t)=0&t\gt0\\ z(1,t)=0&t\gt0 \end{cases} [/math]

Resolviendo el sistema por separación de variables, obtenemos como solución [math] z(x,t)=\sum_{k=1}^\infty\frac{2(-1)^{k+1}}{k\pi}\sin(k\pi x)e^{-k^2\pi^2 t} [/math]. Deshaciendo los cambios que hemos realizado, obtenemos la solución al problema original, [math] u(x,t)=\sum_{k=1}^\infty(-1)\frac{2(-1)^{k+1}}{k\pi}\sin(k\pi x)e^{-k^2\pi^2 t}+x [/math]. Representamos la solución [math] u(x,t) [/math] para [math]x \in [0,1] [/math] en el intervalo de tiempo [math]t \in [0,1] [/math] tomando solo los 10 primeros términos de la serie.

Representación de los [math] 10 [/math] primeros términos de la solución [math] u(x,t) [/math] para [math] (x,t)\in [0,1] \times [0,1] [/math]

En la imagen, podemos ver como a medida que agregamos términos, la solución se vuelve más precisa. Además, vamos a comprobar como las condiciones iniciales de nuestro sistema concuerdan con la solución obtenida. En la posición [math] x=0 [/math] la solución es la constante [math] 0 [/math] y en la posición [math] x=1[/math] la solución es constantemente [math] 1 [/math] . Sin embargo, a medida que aumenta la posición [math] x [/math] en el tiempo inicial, la solución no vale constantemente 0, sino que oscila temporalmente. Esto es debido a que la solución original no es una función continua cerca de [math](x,t)=(1,0) [/math]. Por un lado, tenemos la condición [math]u(x,0)=0[/math], mientras que por el otro lado [math]u(1,t)=1[/math]. Debido a la incapacidad de las series de Fourier para representar de manera precisa las discontinuidades, las sumas parciales de Fourier tienden a "sobrepasar" la discontinuidad en un intento de ajustarse a la función. Este suceso es conocido como el fenómeno de Gibbs y es el responsable de las fluctuaciones mencionadas.

A continuación, aumentamos el valor de la serie hasta [math]k=1000[/math] para ver cómo varia:

Solución [math]u(x,t) [/math], para [math] t \in [0,1][/math] tomando los k primeros términos de la solución.

Como se puede apreciar en esta imagen, cuando [math] k=1000 [/math], este último caso no se aprecian estas oscilaciones. A medida que se agregan más términos a la serie, las oscilaciones disminuyen en amplitud.

Además, esta discontinuidad en el instante inicial provoca que la función no sea diferenciable en [math] t=0 ~s [/math]. De hecho, si calculamos la derivada espacial en la solución término a término de manera analítica en dicho instante, obtenemos que [math] \sum_{k=1}^{\infty} \frac{d}{dx} u_{k}(x,0) + \frac{d}{dx} x=1+\sum_{k=1}^{\infty} 2(-1)^{k} cos(k\pi x) [/math]. En particular, si evaluamos en [math] x=0 [/math]

[math] 1+\sum_{k=1}^{\infty} \frac{d}{dx} u_{k}(0,0) =1+\sum_{k=1}^{\infty} 2(-1)^{k} [/math]

Esta es una serie que no converge, pues toma los valores [math] -1 [/math] y [math] 1 [/math] alternadamente. De manera equivalente, evaluando en [math] x=1 [/math]

[math] 1+\sum_{k=1}^{\infty} u_{k ~ x}(1,0) =1+\sum_{k=1}^{\infty} 2(-1)^{k} cos(k\pi) =1+\sum_{k=1}^{\infty} 2[/math]

Esta última serie tampoco converge, pues el límite de sus términos es claramente no nulo.

Sin embargo, para instantes de tiempo no nulos la solución converge uniformemente y, por ello, se verifica que

[math] u_x(x,t)=\sum_{k=1}^{\infty} \frac{d}{dx} u_{k}(x,t) + \frac{d}{dx} x, ~\quad t\gt0 [/math]

Por lo tanto, concluimos que no tiene sentido analítico calcular la derivada espacial en [math] t=0 [/math].

Asimismo, se ve con claridad como la solución obtenida se aproxima al estado estacionario, pues converge a la función [math] v(x)=x [/math] cuando el tiempo es suficientemente grande.

3.1 Flujo del calor

Ahora vamos a analizar cómo varía el flujo del calor entrante y saliente en los extremos de la varilla a lo largo del tiempo. Para ello, conocemos que el flujo sigue la ecuación [math] \overrightarrow{q} = - \nabla u(x,t) [/math], teniendo en cuenta que el coeficiente de difusión vale [math] 1 [/math] en este caso. Considerando la solución con [math] 10 [/math] y [math] 1000 [/math] términos y evaluando el flujo en Matlab para ambos extremos, obtenemos las siguientes gráficas:

Representación del flujo en los extremos con respecto al tiempo, tomando los k primeros términos de la serie.

Recordamos que en el extremo izquierdo, la temperatura vale constantemente [math] 0°C [/math], mientras que en el extremo derecho queda fijada a [math] 1°C [/math]. Por ello, a medida que va transcurriendo el tiempo, el calor se difunde desde donde más temperatura hay a donde menos, es decir, del extremo izquierdo de la varilla al derecho.

Como en el instante inicial, toda la varilla se encuentra a temperatura [math] 0°C [/math], excepto el extremo derecho, el flujo del extremo izquierdo para instantes próximos al inicial debe valer [math] 0 ~ \frac{cal}{m^2 s} [/math] aproximadamente pues los puntos próximos a él tienen la misma temperatura. Tal y como se cumple cuando aumentamos la precisión para k=1000. A medida que transcurre el tiempo, el valor absoluto de este flujo irá aumentando pues los puntos próximos al extremo izquierdo irán adquiriendo temperatura. Sin embargo, el flujo en el extremo derecho tiene un comportamiento totalmente contrario. En el instante inicial, hay una discontinuidad abrupta de temperatura en el extremo derecho, lo que provoca un flujo cuyo valor absoluto es realmente elevado. Además, el calor se irá transmitiendo lentamente a lo largo de la varilla, provocando que cada vez el valor absoluto del flujo sea más pequeño.

Por otra parte, cabe destacar que, ambos flujos toman valores negativos, lo cual indica que su sentido es opuesto al propio eje direccional, provocando que sea saliente en el extremo izquierdo y entrante en el derecho.

Cabe destacar que, tal y como hemos mencionado antes, la solución no es derivable con respecto al espacio en el instante inicial, y por ello, no tiene sentido estudiar el flujo en [math] t=0 [/math].

Finalmente, es imprescindible mencionar que esta gráfica es excelente a la hora de estudiar la convergencia de la solución a la solución estacionaria. El instante de tiempo en el que ambos flujos se estabilizan tomando el mismo valor simboliza que, en la varilla, el calor que entra por un extremo es el mismo que el que sale por el otro extremo, lo que obviamente es una consecuencia de haber alcanzado el estado estacionario.

3.2 Código de Matlab

Todas las gráficas que se han representado en esta sección, son resultado del siguiente código:

clear all
close all

% COEFICIENTES DE FOURIER
c_k=@(x,k) (-1).*(((2.*(-1).^(k+1))./(k.*pi)));

% TÉRMINO K-ÉSIMO DE LA SOLUCIÓN DEL PROBLEMA HABIENDO RESTADO LA SOLUCIÓN ESTACIONARIA
w_k=@(x,t,k,D,coef_fourier) coef_fourier.*sin(k.*pi.*x).*exp(-D.*(k.^2.*pi.^2.*t)); 

% DEFINIMOS LA SOLUCIÓN ESTACIONARIA
sol_estacionaria=@(x) x;

% VALORES DE LOS COEFICIENTES
coeficiente_conductividad_termica=1;
calor_especifico=1;

% DATOS PARA LA DISCRETIZACIÓN DEL ESPACIO
extremo_izquierdo=0;
extremo_derecho=1;
paso_espacio=0.01;

punto_espacio_estudio=1/2; % punto del espacio concreto (entre extremo_izquierdo y extremo_derecho) que queremos estudiar

% DATOS PARA LA DISCRETIZACIÓN DEL TIEMPO
tiempo_inicial=0;
tiempo_final=1;
paso_tiempo=0.0001;

% DISCRETIZACIÓN DEL ESPACIO Y DEL TIEMPO
[X,T]=meshgrid(extremo_izquierdo:paso_espacio:extremo_derecho,tiempo_inicial:paso_tiempo:tiempo_final); % matrices de mallado espacio-temporal

w=zeros(size(X)); % inicializamos la variable w (solución del problema habiendo restado la solución estacionaria)

% SOLUCIÓN EN EL INTERVALO DE TIEMPO [0,1] TOMANDO LOS k_end PRIMEROS TÉRMINOS DE LA SERIE 
k_end=1000; % último valor de k del sumatorio, siendo k el índice del sumatorio de la solución

for ks=1:k_end % queremos resolver el problema para distintos valores de k
    w=w+w_k(X,T,ks,coeficiente_conductividad_termica/calor_especifico,c_k(X(1,:),ks));
    u=w+sol_estacionaria(X); % deshacemos el cambio de variable sumando la solución estacionaria
    if ks <= 10 % representamos las gráficas de a lo sumo los 10 primeros términos de la solución
        figure(1)
        sgtitle("Solución en el intervalo de tiempo ["+num2str(tiempo_inicial)+","+num2str(tiempo_final)+"] tomando los k primeros términos de la serie")
        subplot(2,min(ceil(k_end/2),5),ks)
        hold on
        mesh(X,T,u)
        xlabel('Espacio')
        ylabel('Tiempo')
        titulo="k = "+num2str(ks);
        title(titulo)
        colormap('jet')
        colorbar
    end

    % COMPARACIÓN DE LA SOLUCIÓN PARA k=10 Y k=k_end
    if (k_end>10) & (ks==10)
        figure(2)
        sgtitle("Solución en el intervalo de tiempo ["+num2str(tiempo_inicial)+","+num2str(tiempo_final)+"] tomando los k primeros términos de la serie")
        subplot(1,2,1)
        mesh(X,T,u)
        colorbar
        colormap('jet')
        xlabel('Espacio')
        ylabel('Tiempo')
        title("k = "+num2str(ks))

        figure(3) % CALCULAR EL FLUJO ENTRANTE Y SALIENTE EN AMBOS EXTREMOS A LO LARGO DEL TIEMPO
        sgtitle('Representación del flujo con respecto al tiempo en los extremos para distintos valores de k')
        subplot(1,2,1)
        deriv=gradient(u,paso_espacio); % derivada de la solución con respecto al espacio
        flujo0=-deriv(:,1); % calculamos el flujo en el extremo izquierdo
        flujo1=-deriv(:,end);% calculamos el flujo en el extremo derecho
        plot(T(:,1),flujo0,'r','linewidth',1.5)
        hold on
        plot(T(:,1),flujo1,'b','linewidth',1.5)
        xlabel('Tiempo')
        ylabel('Flujo')
        legend('Flujo en el extremo izquierdo','Flujo en el extremo derecho','Location','Southeast')
        title("k="+num2str(ks))


    elseif (ks==k_end)
        figure(2)
        if (ks>10)
            subplot(1,2,2)
        end
        mesh(X,T,u)
        colorbar
        xlabel('Espacio')
        ylabel('Tiempo')
        titulo="k = "+num2str(ks);
        title(titulo)

        figure(3) % CALCULAR EL FLUJO ENTRANTE Y SALIENTE EN AMBOS EXTREMOS A LO LARGO DEL TIEMPO
        if (ks>10)
            subplot(1,2,2)
        end
        deriv=gradient(u,paso_espacio); % derivada de la solución con respecto al espacio
        flujo0=-deriv(:,1); % calculamos el flujo en el extremo izquierdo
        flujo1=-deriv(:,end);% calculamos el flujo en el extremo derecho
        plot(T(:,1),flujo0,'r','linewidth',1.5)
        hold on
        plot(T(:,1),flujo1,'b','linewidth',1.5)
        axis([0,1,-100,0])
        xlabel('Tiempo')
        ylabel('Flujo')
        legend('Flujo en el extremo izquierdo','Flujo en el extremo derecho','Location','Southeast')
        title("k="+num2str(ks))
    end
end

En el código, lo primero que hemos hecho ha sido definir los coeficientes de Fourier, c_k, y el término k-ésimo de la solución del problema después de haber restado la solución estacionaria,w_k. Además, se define la solución estacionaria del problema, sol_estacionaria.

Luego hemos establecido valores para los coeficientes de conductividad térmica y el calor específico y se ha discretizado tanto el tiempo como el espacio para crear una malla espaciotemporal X y T mediante la función meshgrid. De esta manera, pasamos de un problema continuo a un problema discreto con un conjunto finito de puntos. Esta discretización nos permite representar la función gráficamente y utilizar métodos numéricos para el cálculo de los coeficientes de Fourier, como veremos posteriormente.

A continuación, calculamos la solución w homogeneizada con k-ésimos términos [math] (k=1,\dots,100) [/math], mostrando en gráficas las distintas soluciones y su evolución con respecto al espacio y tiempo. Además, graficamos la solución para [math] k=10 [/math] y [math] k=1000 [/math] y calculamos los flujos en los extremos, empleando la función gradient. Aunque analíticamente no es posible calcular dichos flujos en el instante inicial, hemos hecho uso de esta función para calcular la derivada de manera numérica. Esto permite que el lector no se vea obligado a calcular dicha derivada por sí mismo, pero puede cometer errores como esta aproximación, así siempre debemos corroborar los resultados obtenidos en estas gráficas con la base teórica.

Cabe destacar que todos estos datos están al comienzo del programa con la idea de que sean fácilmente identificables y modificables para el lector, y así poder estudiar distintas situaciones iniciales que no abordamos en este trabajo.

3.3 Simulación

Para comprender mejor el comportamiento de la varilla, hemos realizado una simulación en la que se representa una varilla de longitud [math] 1 m [/math] y se puede apreciar cómo se transmite el calor a lo largo del tiempo. Para ello, hemos elaborado el siguiente código de Matlab:

clear all
close all
 
% Coeficientes de Fourier
c_k=@(x,k) (-1).*(((2.*(-1).^(k+1))./(k.*pi)));
 
% Término k-ésimo de la solución del problema habiendo restado la solución estacionaria
w_k=@(x,t,k,D,coef_fourier) coef_fourier.*sin(k.*pi.*x).*exp(-D.*(k.^2.*pi.^2.*t)); 
 
% Definimos la solución estacionaria
sol_estacionaria=@(x) x;
 
% Valores de los coeficientes
coeficiente_conductividad_termica=1;
calor_especifico=1;
 
% Definir la longitud y radio de la varilla
longitud = 1;
radio = 0.25;
 
% Definir el número de puntos en la varilla
num_puntos_longitudinal = 1001; % Número de puntos a lo largo de la longitud
num_puntos_circunferencia = 200; % Número de puntos alrededor de la circunferencia
 
% Definir el tiempo y el paso de tiempo
tiempo_final = 1; 
paso_tiempo = 0.005;
tiempo = 0:paso_tiempo:tiempo_final;
 
% Generar la malla de puntos en la varilla
longitudinal = linspace(0, longitud, num_puntos_longitudinal);
circunferencia = linspace(0, 2*pi, num_puntos_circunferencia);
[Longitudinal, Circunferencia] = meshgrid(longitudinal, circunferencia);
 
% Crear un vídeo
pelicula=VideoWriter('Ej1_Ap1_5_Tuberia','MPEG-4'); %creo el video
pelicula.FrameRate=14; %controla velocidad
open(pelicula); %abrimos el vídeo
 
% Bucle para cada paso de tiempo
for t = tiempo
title(['Varilla en t = ' num2str(t)]);
    figura=figure(1);
    
    % Calcular la temperatura en cada punto de la varilla en el tiempo t
    w = zeros(size(Longitudinal));
    for k = 1:1000
        w = w + w_k(Longitudinal,t,k,coeficiente_conductividad_termica/calor_especifico,c_k(Longitudinal(1,:),k));
    end
    u=w+sol_estacionaria(Longitudinal);
 
    % Graficar la varilla con los colores según la temperatura
    colormap(jet); % Mapa de colores de azul (frío) a rojo (caliente)
    surf(Longitudinal, radio * cos(Circunferencia), radio * sin(Circunferencia), u,'EdgeColor','none');
    colorbar;
    clim([0, 1]); % Fijar el rango de colores
    
    % Añadir etiquetas y ajustes de la figura
    title(['Varilla en t = ' num2str(t)]);
    xlabel('Longitud');
    axis([0 longitud -radio radio -radio radio]);
    view(3);
    
    % Añadir el frame al vídeo
    imagen=getframe(figura);
    writeVideo(pelicula,imagen); %inserto la imagen
    
end
 
% Cerrar el vídeo
close(pelicula);

Al igual que en el otro código elaborado en esta sección, lo primero que definimos es la solución estacionaria del problema, así como los coeficientes de Fourier y los términos de la solución homogeneizada. Además, lo siguiente que debemos determinar son el coeficiente de conductividad térmica y el calor específico. En este caso, tomaremos ambos valores iguales a [math] 1 [/math]. Sin embargo, en la próxima sección estudiaremos el efecto que tienen estos valores en la solución.

A continuación, pasamos a introducir los datos que nos ayudarán a la representación de la varilla metálica, como son la longitud y el radio de esta o el número de puntos que pondremos a lo largo del eje x y el número de puntos que tomaremos para formar las circunferencias alrededor del eje longitudinal. Con estos datos generamos la malla de puntos de la varilla. Por otra parte, es fundamental definir los tiempos inicial y final, así como la longitud del intervalo temporal.

Para poder almacenar el resultado en un vídeo, hacemos uso de los comandos correspondientes propios de Matlab, entre los que destacan VideoWriter y FrameRate. Para cada uno de los tiempos definidos, debemos sumar los k-términos de la solución, siendo k el número que deseemos. A este resultado debemos añadirle la solución estacionaria para deshacer el cambio de variable y almacenamos el resultado en la variable u, que representa la temperatura de cada punto de la malla. Esta variable nos permitirá definir la escala de colores con la que vamos a representar la varilla metálica, por lo que a continuación, representamos la varilla haciendo uso de los datos de temperatura en cada punto para asignarle un color. Finalmente, añadimos la imagen generada al vídeo que se está creando y repetimos el proceso hasta haber recorrido todos los instantes de tiempo.

Cabe destacar que, el código realizado puede dar error en su ejecución si no se cuenta con una versión de Matlab de [math] 2022 [/math] o superior. En caso de no contar con dicha versión, basta sustituir en la línea 43 del código clim por caxis.

El resultado obtenido con este código es la siguiente simulación:

Simulación que representa cómo se distribuye el calor a lo largo de la varilla a medida que transcurre el tiempo

En esta simulación, apreciamos como la temperatura, partiendo de [math]0 °C [/math], aumenta de derecha a izquierda de acuerdo con las condiciones impuestas en el sistema de EDP’s y acorde a lo deducido del comportamiento del flujo en las gráficas anteriores.

4 Cambio de coeficientes

Planteamos ahora el mismo escenario anterior pero con coeficiente de conductividad térmica distintos: [math] k_1 =\frac{1}{2} \frac{cal}{ºC s m}, k_2=\frac{1}{6} \frac{cal}{ºC s m}[/math], en lugar de [math] k=1 \frac{cal}{ºC s m} [/math]. El calor específico [math] c [/math] continuará siendo igual a [math] 1 \frac{cal}{ºC g} [/math]. Por tanto, los sistemas quedarían como:

[math] \left \{ \begin{array}{ll} \frac{\partial u}{\partial t}-\frac{ki}{c} \frac{\partial ^2 u}{\partial x^2}=0 & \quad 0 \lt x \lt 1, 0 \lt t \lt 1, i = 1,2 \\ u(x,0)=0, & \quad 0 \lt x \lt 1, \\ u(0,t)=0, & \quad 0 \lt t \lt 1, \\ u(1,t)=1 & \quad 0 \lt t \lt 1. \end{array} \right. [/math]

De manera análoga a los casos anteriores, obtenemos que la solución es de la forma

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

En primer lugar, es imprescindible recordar para la comparación que cuanto mayor sea el coeficiente de conductividad térmica [math] k [/math], mayor será la capacidad del material para conducir el calor. Además, estudiaremos cómo se comporta la solución del sistema homogeneizado en el punto intermedio de la barra, en [math] x = \frac{1}{2} [/math], y la compararemos para los coeficientes de conductividad térmica [math] k=1 [/math] y [math] k_1=\frac{1}{2} [/math]. Notar que, la solución del sistema homogeneizado equivale a [math] u \left(\frac{1}{2}, t\right) –v\left(\frac{1}{2}\right) [/math].

Para ello, hemos realizado el código presentado a continuación.

clear all
close all
 
% COEFICIENTES DE FOURIER
c_k=@(x,k) (-1).*(((2.*(-1).^(k+1))./(k.*pi)));
 
% TÉRMINO K-ÉSIMO DE LA SOLUCIÓN DEL PROBLEMA HABIENDO RESTADO LA SOLUCIÓN ESTACIONARIA
w_k=@(x,t,k,D,coef_fourier) coef_fourier.*sin(k.*pi.*x).*exp(-D.*(k.^2.*pi.^2.*t)); 
 
% DEFINIMOS LA SOLUCIÓN ESTACIONARIA
sol_estacionaria=@(x) x;
 
% VALORES DE LOS COEFICIENTES
coeficiente_conductividad_termica=[1,1/2,1/6];
calor_especifico=[1,1,1];
 
% DATOS PARA LA DISCRETIZACIÓN DEL ESPACIO
extremo_izquierdo=0;
extremo_derecho=1;
paso_espacio=0.01;
 
punto_espacio_estudio=1/2; % punto del espacio concreto (entre extremo_izquierdo y extremo_derecho) que queremos estudiar
 
% DATOS PARA LA DISCRETIZACIÓN DEL TIEMPO
tiempo_inicial=0;
tiempo_final=1;
paso_tiempo=0.0001;
 
% DISCRETIZACIÓN DEL ESPACIO Y DEL TIEMPO
[X1,T]=meshgrid(extremo_izquierdo:paso_espacio:extremo_derecho,tiempo_inicial:paso_tiempo:tiempo_final); % matrices de mallado espacio-temporal
X1(:,end+1)=punto_espacio_estudio*ones(length(X1(:,1)),1); % vector con el punto del espacio a estudiar
X=transpose(sortrows(transpose(X1))); % añadimos el punto del espacio a estudiar en la matriz de mallado espacial
T(:,end+1)=T(:,end); % redimensionamos la matriz de mallado temporal
 
% SOLUCIÓN EN EL INTERVALO DE TIEMPO [0,1] TOMANDO LOS k_end PRIMEROS TÉRMINOS DE LA SERIE 
k_end=200; % último valor de k del sumatorio, siendo k el índice del sumatorio de la solución
 
for coefi=1:length(coeficiente_conductividad_termica)
    w_coef=zeros(size(X)); % inicializamos la variable w_coef1 (solución del problema con el primer coeficiente habiendo restado la solución estacionaria)
    for ks=1:k_end % resolvemos el problema con los primeros k_end términos.
        w_coef=w_coef+w_k(X,T,ks,coeficiente_conductividad_termica(coefi)/calor_especifico(coefi),c_k(X(1,:),ks)); 
    end
    u=w_coef+sol_estacionaria(X); % deshacemos el cambio de variable sumando la solución estacionaria
    
    % REPRESENTACIONES DE LAS SOLUCIONES PARA LOS DOS COEFICIENTES DE CONDUCTIVIDAD
    figure(1)
    subplot(1,length(coeficiente_conductividad_termica),coefi)
    mesh(X,T,u)
    colormap('jet')
    colorbar
    xlabel('Espacio')
    ylabel('Tiempo')
    title("D = "+num2str(coeficiente_conductividad_termica(coefi)/calor_especifico(coefi)))
    sgtitle("Representación de los primeros "+num2str(k_end)+" términos de la solución según los valores de D")
 
    % REPRESENTAMOS LA DIFERENCIA ENTRE LA SOLUCIÓN Y LA SOLUCIÓN ESTACIONARIA EN punto_espacio_estudio
    figure(2)
    plot(T(:,1),w_coef(:,ceil(length(w_coef(1,:))/2)),'linewidth',1.5,'DisplayName',"Solución con D="+num2str(coeficiente_conductividad_termica(coefi)/calor_especifico(coefi)))
    hold on
    xlabel('Tiempo')
    ylabel('Temperatura')
    title("Diferencia entre la solución y la solución estacionaria en x="+num2str(punto_espacio_estudio)+" con respecto al tiempo")
 
end
 
figure(2)
legend('Location','southeast')

De nuevo, en primer lugar se definen los coeficientes de Fourier de la extensión impar de la función [math] x [/math] (solución estacionaria) como [math] c_k [/math], además del k-ésimo término de la solución homogeneizada y la solución estacionaria del sistema.

Posteriormente se incluyen distintos datos o parámetros modificables según las condiciones que deseemos establecer. Entre ellos se encuentran los valores de los coeficientes de conductividad térmica [math] k [/math], en forma de vector, que busquemos comparar.

Por otro lado, tenemos el dominio espaciotemporal que busquemos establecer. Al igual que en todos los códigos anteriores, de cara a aproximar los resultados continuos de una manera discreta, debemos establecer una malla espaciotemporal. Cuantos más puntos tenga esta malla, mayor precisión tendrá la aproximación al caso continuo.

Es destacable la definición de punto_espacio_estudio. Este punto es un punto específico dentro del dominio espacial en el cual queremos realizar el estudio. En dicho punto calcularemos cómo evoluciona la solución homogeneizada a medida que avanza el tiempo. Para añadir este punto estudio a la malla espaciotemporal, añadimos una columna con el punto repetido en todas las filas a dicha matriz de mallado, y posteriormente ordenamos la matriz por columnas.

Por último, realizamos una construcción análoga a la ya mencionada en el código de la sección anterior, a diferencia de que el proceso se debe repetir para cada uno de los coeficientes [math] k [/math] para los que deseemos conocer la solución.

Representación de los [math] 200 [/math] primeros términos de la solución para los coeficientes de conductividad térmica [math] k=1 \frac{cal}{ºC s m} [/math] y [math] k_1=\frac{1}{2} \frac{cal}{ºC s m} , k_2=\frac{1}{6} \frac{cal}{ºC s m}[/math] en el intervalo de tiempo [math] [0,1] [/math].

Con respecto a esta primera imagen, como se puede observar, los dos primeros casos son muy similares, aunque hayamos cambiado el coeficiente de conductividad térmica ligeramente. Sin embargo, si comparamos los casos [math] k=1 \frac{cal}{ºC s m} [/math] y [math] k_2=\frac{1}{6} \frac{cal}{ºC s m}[/math], sí se observa una diferencia en el comportamiento. La superficie representada con [math] k_2= \frac{1}{6} \frac{cal}{ºC s m} [/math] presenta una mayor curvatura en el instante de tiempo final representado ([math] t=1 s[/math]).

Esto se debe a que aunque las tres soluciones se aproximan a la solución estacionaria [math] v(x,t)=x [/math], el coeficiente de conductividad térmica influye en la velocidad de convergencia de la solución hacia la solución estacionaria. Por ello, para un instante de tiempo [math] t=1 s[/math] se observan diferentes curvaturas en las superficies. En caso de aumentar el intervalo de tiempo lo suficiente, las diferencias en el instante final serían despreciables. Se anima al lector a cambiar el valor del parámetro tiempo_final por [math] 3’5 s[/math] en el código recién explicado.

Representación de la solución del sistema homogeneizado en [math] x=\frac{1}{2} [/math] con los coeficientes de conductividad térmica [math] k=1 \frac{cal}{ºC s m}[/math] en color azul, [math] k_1=\frac{1}{2} \frac{cal}{ºC s m}[/math] en color rojo y [math] k_2=\frac{1}{6} \frac{cal}{ºC s m}[/math] en color amarillo.

Esta conclusión se puede observar perfectamente en esta segunda gráfica, en la que se muestra la solución del sistema homogeneizado en [math] x=\frac{1}{2} [/math] con los diferentes coeficientes de conductividad térmica. Como en la gráfica se representa la diferencia entre la solución y la estacionaria, cuanto más cercanos a cero sean estos valores, menos tiempo quedará para la convergencia de la solución.

Si prestamos atención a la gráfica, se puede observar a simple vista como, a medida que disminuimos los valores de los coeficientes de conductividad térmica [math] k [/math], el tiempo necesario para llegar a una temperatura nula es cada vez mayor. Luego, la velocidad de convergencia de la solución a la estacionaria aumenta a medida que aumenta el coeficiente de conductividad térmica.

Todo ello tiene sentido con el concepto de coeficiente de conductividad térmica dado que, los coeficientes de conductividad térmica altos conducen mejor el calor y, por tanto, tardan una menor cantidad de tiempo en estabilizar la temperatura del cuerpo con respecto a la temperatura ambiental.

4.1 Simulación

En esta sección, también se ha elaborado un código para la simulación y visualización del problema en el caso de la varilla metálica.

clear all
close all
 
% Coeficientes de Fourier
c_k=@(x,k) (-1).*(((2.*(-1).^(k+1))./(k.*pi)));
 
% Término k-ésimo de la solución del problema habiendo restado la solución estacionaria
w_k=@(x,t,k,D,coef_fourier) coef_fourier.*sin(k.*pi.*x).*exp(-D.*(k.^2.*pi.^2.*t)); 
 
% Definimos la solución estacionaria
sol_estacionaria=@(x) x;
 
% Valores de los coeficientes
coeficiente_conductividad_termica=[1,1/2,1/6];
calor_especifico=[1,1,1];
 
% Definir la longitud y radio de la varilla
longitud = 1;
radio = 0.25;
 
% Definir el número de puntos en la varilla
num_puntos_longitudinal = 1001; % Número de puntos a lo largo de la longitud
num_puntos_circunferencia = 200; % Número de puntos alrededor de la circunferencia
 
% Definir el tiempo y el paso de tiempo
tiempo_final = 1; 
paso_tiempo = 0.005;
tiempo = 0:paso_tiempo:tiempo_final;
 
% Generar la malla de puntos en la varilla
longitudinal = linspace(0, longitud, num_puntos_longitudinal);
circunferencia = linspace(0, 2*pi, num_puntos_circunferencia);
[Longitudinal, Circunferencia] = meshgrid(longitudinal, circunferencia);
 
% Crear un vídeo
pelicula=VideoWriter('Ej1_Ap6_Tuberia','MPEG-4'); %creo el video
pelicula.FrameRate=14; %controla velocidad
open(pelicula); %abrimos el vídeo
 
% Bucle para cada paso de tiempo
for t = tiempo
    figura=figure(1);
    sgtitle("Varilla en t = "+num2str(t)+" para distintos valores de D")
    figura.Position=[100,100,1550,550]
    for coefi=1:length(coeficiente_conductividad_termica)
    
        % Calcular la temperatura en cada punto de la varilla en el tiempo t
        w = zeros(size(Longitudinal));
        for k = 1:1000
            w = w + w_k(Longitudinal,t,k,coeficiente_conductividad_termica(coefi)/calor_especifico(coefi),c_k(Longitudinal(1,:),k));
        end
        u=w+sol_estacionaria(Longitudinal);
    
        % Graficar la varilla con los colores según la temperatura
        subplot(1,length(coeficiente_conductividad_termica),coefi)
        colormap(jet); % Mapa de colores de azul (frío) a rojo (caliente)
        surf(Longitudinal, radio * cos(Circunferencia), radio * sin(Circunferencia), u,'EdgeColor','none');
        colorbar;
        clim([0, 1]); % Fijar el rango de colores
        
        % Añadir etiquetas y ajustes de la figura
        title("D="+coeficiente_conductividad_termica(coefi)/calor_especifico(coefi));
        xlabel('Longitud');
        axis([0 longitud -radio radio -radio radio]);
        view(3);
        
        
    end
    % Añadir el frame al vídeo
    imagen=getframe(figura);
    writeVideo(pelicula,imagen); %inserto la imagen
end
 
% Cerrar el vídeo
close(pelicula);

La construcción del código es análoga a la ya mencionada en la sección 3.3, pero con una pequeña diferencia. Se define un vector de coeficientes de conductividad térmica [math] k [/math] y se repite el proceso de representación para cada uno de los valores.

Simulación de la solución en la varilla metálica. Nota: Si no se visualiza la animación de este archivo gif, pinchar en el archivo para su visualización.

En este caso, se aprecian perfectamente las observaciones realizadas. El coeficiente de conductividad térmica afecta principalmente a la velocidad con la que se transmite el calor. Si observamos el primer y tercer video, la velocidad a la que se transmite el calor es muy distinta. En el primer caso, con coeficiente de difusión [math] D=1 [/math], el sistema alcanza el estado estacionario mucho más rápido. De hecho, el caso en el que el coeficiente de difusión [math] D=1/6 [/math], el sistema aún no ha alcanzado el estado estacionario cuando finaliza el vídeo.

5 Cambio de condiciones

Una vez hemos estudiado cómo influye el realizar un cambio en el coeficiente de conductividad térmica en el problema original, nos planteamos cómo varía la solución si suponemos que la temperatura en el extremo derecho es también de [math] 0 °C[/math] pero inicialmente la temperatura de la varilla viene dada por la función [math] u_0(x)=max\{0,1-4|x-1/2|\}[/math].

Por tanto, tendríamos el siguiente sistema:

[math] \left \{ \begin{array}{ll} \frac{\partial u}{\partial t}- \frac{\partial ^2 u}{\partial x^2}=0 & \quad 0 \lt x \lt 1, 0 \lt t \lt 1\\ u(x,0)= u_0(x)=max\{0,1-4|x-1/2|\}, & \quad 0 \lt x \lt 1, \\ u(0,t)=0, & \quad 0 \lt t \lt 1, \\ u(1,t)=0 & \quad 0 \lt t \lt 1. \end{array} \right. [/math]

La solución general de este sistema permanece igual que en sistema original a diferencia de que, en este caso, al tener las condiciones frontera homogeneizadas, la solución particular es nula. Es decir, la solución general es de la forma:

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

La única diferencia es la determinación de los coeficientes de Fourier [math] c_k [/math], en este caso, resulta mucho más sencillo aproximar el valor de las integrales de dichos coeficientes por la fórmula del trapecio.

A continuación, se presenta un código en Matlab que sigue la misma estructura que el presentado en las secciones [math] 3.2 [/math]. A diferencia de este, al haber cambiado la condición inicial por la función [math] u_0(x) [/math] y la condición frontera en el extremo derecho a una temperatura inicial nula, únicamente cambia que sol_estacionaria(x) es un vector de ceros y la determinación de los coeficientes de Fourier, en este caso calculados por la fórmula del trapecio.

clear all
close all
 
% COEFICIENTES DE FOURIER CALCULADOS MEDIANTE LA APROXIMACIÓN POR LA FÓMULA DEL TRAPECIO
c_k=@(x,k) trapz(x,max(0,1-4.*abs(x-1/2*ones(size(x)))).*sin(k*pi.*x))./trapz(x,(sin(k*pi.*x)).^2);
 
% TÉRMINO K-ÉSIMO DE LA SOLUCIÓN DEL PROBLEMA HABIENDO RESTADO LA SOLUCIÓN ESTACIONARIA
w_k=@(x,t,k,D,coef_fourier) coef_fourier.*sin(k.*pi.*x).*exp(-D.*(k.^2.*pi.^2.*t)); 
 
% DEFINIMOS LA SOLUCIÓN ESTACIONARIA
sol_estacionaria=@(x) zeros(size(x));
 
% VALORES DE LOS COEFICIENTES
coeficiente_conductividad_termica=1;
calor_especifico=1;
 
% DATOS PARA LA DISCRETIZACIÓN DEL ESPACIO
extremo_izquierdo=0;
extremo_derecho=1;
paso_espacio=0.001;
 
punto_espacio_estudio=1/2; % punto del espacio concreto (entre extremo_izquierdo y extremo_derecho) que queremos estudiar
 
% DATOS PARA LA DISCRETIZACIÓN DEL TIEMPO
tiempo_inicial=0;
tiempo_final=1;
paso_tiempo=0.001;
 
% DISCRETIZACIÓN DEL ESPACIO Y DEL TIEMPO
[X,T]=meshgrid(extremo_izquierdo:paso_espacio:extremo_derecho,tiempo_inicial:paso_tiempo:tiempo_final); % matrices de mallado espacio-temporal
 
w=zeros(size(X)); % inicializamos la variable w (solución del problema habiendo restado la solución estacionaria)
 
% SOLUCIÓN EN EL INTERVALO DE TIEMPO [0,1] TOMANDO LOS k_end PRIMEROS TÉRMINOS DE LA SERIE 
k_end=400; %último valor de k del sumatorio, siendo k el índice del sumatorio de la solución
 
for ks=1:k_end % queremos resolver el problema para distintos valores de k
    w=w+w_k(X,T,ks,coeficiente_conductividad_termica/calor_especifico,c_k(X(1,:),ks));
    u=w+sol_estacionaria(X); % deshacemos el cambio de variable sumando la solución estacionaria
    if ks <= 10 % representamos las gráficas de a lo sumo los 10 primeros términos de la solución
        figure(1)
        subplot(2,min(ceil(k_end/2),5),ks)
        hold on
        mesh(X,T,u)
        colormap('jet')
        colorbar
        xlabel('Espacio')
        ylabel('Tiempo')
        titulo="k = "+num2str(ks);
        title(titulo)
    end
    
    % COMPARACIÓN DE LA SOLUCIÓN PARA k=10 Y k=k_end
    if (k_end>10) & (ks==10)
        figure(2)
        sgtitle("Solución en el intervalo de tiempo ["+num2str(tiempo_inicial)+","+num2str(tiempo_final)+"] tomando los k primeros términos de la serie")
        subplot(1,2,1)
        mesh(X,T,u)
        colormap('jet')
        colorbar
        xlabel('Espacio')
        ylabel('Tiempo')
        titulo="k = "+num2str(ks);
        title(titulo)
 
        figure(3) % CALCULAR EL FLUJO ENTRANTE Y SALIENTE EN AMBOS EXTREMOS A LO LARGO DEL TIEMPO
        sgtitle('Representación del flujo con respecto al tiempo en los extremos para distintos valores de k')
        subplot(1,2,1)
        deriv=gradient(u,paso_espacio); % derivada de la solución con respecto al espacio
        flujo0=-deriv(:,1); % calculamos el flujo en el extremo izquierdo
        flujo1=-deriv(:,end);% calculamos el flujo en el extremo derecho
        plot(T(:,1),flujo0,'r','linewidth',1.5)
        hold on
        plot(T(:,1),flujo1,'b','linewidth',1.5)
        xlabel('Tiempo')
        ylabel('Flujo')
        legend('Flujo en el extremo izquierdo','Flujo en el extremo derecho','Location','Southeast')
        title("k="+num2str(ks))
 
 
    elseif (ks==k_end)
        figure(2)
        if (ks>10)
            subplot(1,2,2)
        end
        mesh(X,T,u)
        colormap('jet')
        colorbar
        xlabel('Espacio')
        ylabel('Tiempo')
        titulo="k = "+num2str(ks);
        title(titulo)
 
        figure(3) % CALCULAR EL FLUJO ENTRANTE Y SALIENTE EN AMBOS EXTREMOS A LO LARGO DEL TIEMPO
        if (ks>10)
            subplot(1,2,2)
        end
        deriv=gradient(u,paso_espacio); % derivada de la solución con respecto al espacio
        flujo0=-deriv(:,1); % calculamos el flujo en el extremo izquierdo
        flujo1=-deriv(:,end);% calculamos el flujo en el extremo derecho
        plot(T(:,1),flujo0,'r','linewidth',1.5)
        hold on
        plot(T(:,1),flujo1,'b','linewidth',1.5)
        xlabel('Tiempo')
        ylabel('Flujo')
        legend('Flujo en el extremo izquierdo','Flujo en el extremo derecho','Location','Southeast')
        title("k="+num2str(ks))
    end
end
figure(4)
scatter(X(1,:),u(1,:),5, u(1,:),'filled')
colormap('jet')
title("Solución para k="+ num2str(k_end)+ " y t=0")
xlabel('Espacio')
ylabel('Temperatura')

Representación de los [math] 10 [/math] primeros términos de la solución con condiciones de frontera nulas y condición inicial la función [math] u_0(x)=max\{0,1-4|x-1/2|\} [/math]representada en el intervalo de tiempo [math] [0,1] [/math].

Cómo se puede observar, la solución obtenida en la primera representación cumple con las condiciones establecidas en el sistema. Además, a medida que aumenta la precisión de la solución (introduciendo un mayor número de términos en la suma), esta alcanza un máximo de temperatura mayor y en un intervalo espacial más pequeño con respecto al punto medio de la varilla. Por otro lado, es destacable como la solución es estrictamente positiva si [math] t\gt0 [/math]. Esto se puede interpretar como que la difusión transporta energía con velocidad infinita de propagación.

Representación de la solución para los k=10 primeros términos (primera imagen) y para los primeros k=400 (segunda imagen)

La mejor manera de entender este fenómeno es observando esta segunda gráfica, en la que se ha representado la solución para los [math] k=10 [/math] y los [math] k=400 [/math] primeros términos. La solución para [math] k=400 [/math] tiene forma de pico. Esto es llamativo dado que la solución no es diferenciable en el instante inicial [math] t=0 [/math].

Representación de la solución para los k=400 primeros términos en el instante de tiempo inicial [math] t=0 [/math]

Sin embargo, en cuanto [math] t [/math] alcanza valores positivos, el calor que se encontraba acumulado en la zona central de la varilla se difunde a los extremos de la misma, provocando que la temperatura en los puntos en los que era nula en el instante inicial, pase a no serla. Cómo esto sucede para cualquier incremento de tiempo tan pequeño como deseemos, esto indica que la velocidad de propagación de calor es infinita.

En relación con la tercera gráfica obtenida con el código citado anteriormente, se puede observar como el flujo en ambos extremos comienza en 0. Esto se debe a que la condición inicial verifica que en los puntos próximos a dichos extremos la temperatura es también cero. Sin embargo, a medida que el calor concentrado alrededor del punto medio se dispersa hacia los extremos de la varilla, el flujo en ambos extremos se ve aumentado en valor absoluto. Cabe destacar que el flujo en el extremo izquierdo de la varilla (representado en color rojo) es negativo y, por ello, saliente; mientras que el flujo en el extremo derecho (en color azul) es positivo y, por ello, también saliente. El hecho de que ambos flujos en los extremos sean salientes indica que la varilla se enfría en su interior. Además, ambos flujos tienden a cero cuando se alcanza la solución estacionaria, lo que indica que la temperatura en la varilla es constante a lo largo de toda ella.

Representación del flujo con respecto al tiempo en los extremos. En color rojo el flujo en el lado izquierdo de la varilla, y en color azul aparece representado el flujo en el lado derecho de la varilla

De nuevo, hemos implementado un código en Matlab de la varilla metálica para mejorar su comprensión. El código realizado es idéntico al realizado en la sección 3.3, únicamente se diferencian en las funciones y condiciones iniciales que las establecemos acordes a este nuevo problema. El código sería el siguiente:

clear all
close all
 
% COEFICIENTES DE FOURIER CALCULADOS MEDIANTE LA APROXIMACIÓN POR LA FÓMULA DEL TRAPECIO
c_k=@(x,k) trapz(x,max(0,1-4.*abs(x-1/2*ones(size(x)))).*sin(k*pi.*x))./trapz(x,(sin(k*pi.*x)).^2);
 
% TÉRMINO K-ÉSIMO DE LA SOLUCIÓN DEL PROBLEMA HABIENDO RESTADO LA SOLUCIÓN ESTACIONARIA
w_k=@(x,t,k,D,coef_fourier) coef_fourier.*sin(k.*pi.*x).*exp(-D.*(k.^2.*pi.^2.*t)); 
 
% DEFINIMOS LA SOLUCIÓN ESTACIONARIA
sol_estacionaria=@(x) zeros(size(x));
 
% VALORES DE LOS COEFICIENTES
coeficiente_conductividad_termica=1;
calor_especifico=1;
 
% Definir la longitud y radio de la varilla
longitud = 1;
radio = 0.25;
 
% Definir el número de puntos en la varilla
num_puntos_longitudinal = 1001; % Número de puntos a lo largo de la longitud
num_puntos_circunferencia = 200; % Número de puntos alrededor de la circunferencia
 
% Definir el tiempo y el paso de tiempo
tiempo_final = 0.55; 
paso_tiempo = 0.001;
tiempo = 0:paso_tiempo:tiempo_final;
 
% Generar la malla de puntos en la varilla
longitudinal = linspace(0, longitud, num_puntos_longitudinal);
circunferencia = linspace(0, 2*pi, num_puntos_circunferencia);
[Longitudinal, Circunferencia] = meshgrid(longitudinal, circunferencia);
 
% Crear un vídeo
pelicula=VideoWriter('Ej1_Ap7_Tuberia','MPEG-4'); %creo el video
pelicula.FrameRate=37; %controla velocidad
open(pelicula); %abrimos el vídeo
 
% Bucle para cada paso de tiempo
for t = tiempo
    figura=figure(1);
    
    % Calcular la temperatura en cada punto de la varilla en el tiempo t
    w = zeros(size(Longitudinal));
    for k = 1:200
        w = w + w_k(Longitudinal,t,k,coeficiente_conductividad_termica/calor_especifico,c_k(Longitudinal(1,:),k));
    end
    u=w+sol_estacionaria(Longitudinal);
 
    % Graficar la varilla con los colores según la temperatura
    colormap(jet); % Mapa de colores de azul (frío) a rojo (caliente)
    surf(Longitudinal, radio * cos(Circunferencia), radio * sin(Circunferencia), u,'EdgeColor','none');
    colorbar;
    clim([0, 1]); % Fijar el rango de colores
    
    % Añadir etiquetas y ajustes de la figura
    title(['Varilla en t = ' num2str(t)]);
    xlabel('Longitud');
    axis([0 longitud -radio radio -radio radio]);
    view(3);
    
    % Añadir el frame al vídeo
    imagen=getframe(figura);
    writeVideo(pelicula,imagen); %inserto la imagen
    if t==0
        imwrite(imagen.cdata,"Ej1_Ap7_Tuberia_Inicio.png") % guardamos la imagen en el instante inicial
    end
    
end
 
% Cerrar el vídeo
close(pelicula);

Representación de la simulación de la varilla metálica para la solución con condiciones de frontera nulas y condición inicial la función [math] u_0(x)=max \{0,1-4|x-1/2| \} [/math] representada en el intervalo de tiempo [math] [0,55] [/math].

Como se puede observar, la varilla comienza con una temperatura inicial más alta en el interior, pero con el trascurso del tiempo, se va enfriando hasta alcanzar el estado estacionario con temperatura nula en toda la varilla.

Representación de la simulación de la varilla metálica para la solución en el instante [math] t=0 [/math]

Por último, dado que la condición inicial establecida (a la que hemos denotado como [math] u_0(x) [/math]), puede resultar algo difícil de visualizar, hemos realizado una representación gráfica de dicha varilla en el instante inicial. Con esta representación gráfica podemos hacernos a la idea de cómo se distribuye la temperatura en [math] t=0 [/math] de una manera más sencilla para así acabar de entender este caso.

6 Cambios de condiciones de frontera

Ahora, en vez de suponer que la temperatura en el extremo derecho es de [math]0ºC [/math], suponemos que la barra está aislada térmicamente en dicho extremo. El sistema de EDP’s quedaría de la siguiente manera:

[math] \begin{cases} u_t-Du_{xx}=0 & x\in(0,1),t\gt0\\ u(x,0)= max\{0,1-4|x-1/2|\}, & 0 \lt x \lt 1, \\ u(0,t)=0&t\gt0\\ u_x(1,t)=0 &t\gt0 \end{cases}[/math]

Resolviendo el sistema por separación de variables obtenemos como solución la serie [math] u(x,t)=\sum_{k=1}^\infty c_{k} sin\left(\left(k \pi-\frac{\pi}{2}\right)x\right)e^{(k \pi -\frac{\pi}{2})^2 Dt }+x [/math]. Esta solución no queda en función de la base trigonométrica usual, no obstante el conjunto de funciones [math]\left\{ sin\left(\left(k \pi-\frac{\pi}{2}\right)x\right)\right\} _{k \in \mathbb{N}} [/math] es un conjunto completo en [math] x\in [0,1] [/math] y es un conjunto ortogonal:

[math] \int^{1}_{0} \sin\left(\left(n \pi-\frac{\pi}{2}\right)x\right) \sin\left(\left(m \pi-\frac{\pi}{2}\right)x\right)dx= \int^{1}_{0} \frac{\cos\left(\left(n\pi-m\pi\right)x\right) -\cos\left(\left(n\pi+m\pi - \pi\right)x\right)}{2}=0[/math], para [math] n,m \in \mathbb{N}[/math]

Por consiguiente, podemos aproximar el dato inicial usando los coeficientes de Fourier asociados a esta nueva base. Calculando estos coeficientes en Matlab mediante la fórmula del trapecio, obtenemos el correspondiente comportamiento de la solución para distinto número de términos de la solución, es decir, valores de k:

clear all
close all

% COEFICIENTES DE FOURIER CALCULADOS MEDIANTE LA APROXIMACIÓN POR LA FÓMULA DEL TRAPECIO
c_k=@(x,k) trapz(x,max(zeros(size(x)),1-4.*abs(x-1/2*ones(size(x)))).*sin((k*pi-pi/2).*x))./trapz(x,(sin((k*pi-pi/2).*x)).^2);

% TÉRMINO K-ÉSIMO DE LA SOLUCIÓN DEL PROBLEMA HABIENDO RESTADO LA SOLUCIÓN ESTACIONARIA
w_k=@(x,t,k,D,coef_fourier) coef_fourier.*sin((k*pi-pi/2).*x).*exp(-D.*((-pi/2+k*pi)^2.*t)); 

% DEFINIMOS LA SOLUCIÓN ESTACIONARIA
sol_estacionaria=@(x) zeros(size(x));

% VALORES DE LOS COEFICIENTES
coeficiente_conductividad_termica=1;
calor_especifico=1;

% DATOS PARA LA DISCRETIZACIÓN DEL ESPACIO
extremo_izquierdo=0;
extremo_derecho=1;
paso_espacio=0.001;

punto_espacio_estudio=1/2; % punto del espacio concreto (entre extremo_izquierdo y extremo_derecho) que queremos estudiar

% DATOS PARA LA DISCRETIZACIÓN DEL TIEMPO
tiempo_inicial=0;
tiempo_final=2;
paso_tiempo=0.001;

% DISCRETIZACIÓN DEL ESPACIO Y DEL TIEMPO
[X,T]=meshgrid(extremo_izquierdo:paso_espacio:extremo_derecho,tiempo_inicial:paso_tiempo:tiempo_final); % matrices de mallado espacio-temporal

w=zeros(size(X)); % inicializamos la variable w (solución del problema habiendo restado la solución estacionaria)

% SOLUCIÓN EN EL INTERVALO DE TIEMPO [0,1] TOMANDO LOS k_end PRIMEROS TÉRMINOS DE LA SERIE 
k_end=1000; %último valor de k del sumatorio, siendo k el índice del sumatorio de la solución

for ks=1:k_end % queremos resolver el problema para distintos valores de k
    w=w+w_k(X,T,ks,coeficiente_conductividad_termica/calor_especifico,c_k(X(1,:),ks));
    u=w+sol_estacionaria(X); % deshacemos el cambio de variable sumando la solución estacionaria
    if ks <= 10 % representamos las gráficas de a lo sumo los 10 primeros términos de la solución
        figure(1)
        subplot(2,min(ceil(k_end/2),5),ks)
        hold on
        mesh(X,T,u)
        colorbar
        colormap('jet')
        xlabel('Espacio')
        ylabel('Tiempo')
        titulo="k = "+num2str(ks);
        title(titulo)
    end

    % COMPARACIÓN DE LA SOLUCIÓN PARA k=10 Y k=k_end
    if (k_end>10) & (ks==10)
        figure(2)
        sgtitle("Solución en el intervalo de tiempo ["+num2str(tiempo_inicial)+","+num2str(tiempo_final)+"] tomando los k primeros términos de la serie")
        subplot(1,2,1)
        mesh(X,T,u)
        colorbar
        colormap('jet')
        xlabel('Espacio')
        ylabel('Tiempo')
        titulo="k = "+num2str(ks);
        title(titulo)

        figure(3) % CALCULAR EL FLUJO ENTRANTE Y SALIENTE EN AMBOS EXTREMOS A LO LARGO DEL TIEMPO
        sgtitle('Representación del flujo con respecto al tiempo en los extremos para distintos valores de k')
        subplot(1,2,1)
        deriv=gradient(u,paso_espacio); % derivada de la solución con respecto al espacio
        flujo0=-deriv(:,1); % calculamos el flujo en el extremo izquierdo
        flujo1=-deriv(:,end);% calculamos el flujo en el extremo derecho
        plot(T(:,1),flujo0,'r','linewidth',1.5)
        hold on
        plot(T(:,1),flujo1,'b','linewidth',1.5)
        xlabel('Tiempo')
        ylabel('Flujo')
        legend('Flujo en el extremo izquierdo','Flujo en el extremo derecho','Location','Southeast')
        title("k="+num2str(ks))

    elseif (ks==k_end)
        figure(2)
        if (ks>10)
            subplot(1,2,2)
        end
        mesh(X,T,u)
        colorbar
        colormap('jet')
        xlabel('Espacio')
        ylabel('Tiempo')
        titulo="k = "+num2str(ks);
        title(titulo)

        figure(3) % CALCULAR EL FLUJO ENTRANTE Y SALIENTE EN AMBOS EXTREMOS A LO LARGO DEL TIEMPO
        if (ks>10)
            subplot(1,2,2)
        end
        deriv=gradient(u,paso_espacio); % derivada de la solución con respecto al espacio
        flujo0=-deriv(:,1); % calculamos el flujo en el extremo izquierdo
        flujo1=-deriv(:,end);% calculamos el flujo en el extremo derecho
        plot(T(:,1),flujo0,'r','linewidth',1.5)
        hold on
        plot(T(:,1),flujo1,'b','linewidth',1.5)
        xlabel('Tiempo')
        ylabel('Flujo')
        legend('Flujo en el extremo izquierdo','Flujo en el extremo derecho','Location','Southeast')
        title("k="+num2str(ks))
    end
end

Representación de la solución para los k primeros término, con k=1,…,10.

En primer lugar, para valores suficientemente altos de k, en el espacio [math]x=0 [/math] la solución se mantiene constantemente nula como indica la condición frontera. Por otra parte, en los espacios intermedios de la varilla para [math] t=0 [/math] la solución toma valores más elevados acorde a la condición inicial [math] u(x,0)= max\{0,1-4|x-1/2|\}, \quad 0 \lt x \lt 1[/math]. Sin embargo, vamos a aumentar el número de términos hasta[math]k=1000[/math] para analizar con mayor precisión el fenómeno que sucede en el extremo derecho

Representación la solución para los 10 y 100 primeros términos.

En comparación con el caso anterior, el calor no se transmite en ambos extremos de la varilla de manera uniforme, sino que hay mayor transporte de calor hacia el lado izquierdo. Esto se debe a que hemos impuesto que el extremo derecho esté aislado térmicamente. Observando la gráfica con detalle, se aprecia que en el extremo izquierdo hay una clara diferencia de temperatura con respecto a los puntos próximos, por lo que en este extremo sí habrá flujo. Sin embargo, cerca del extremo derecho la temperatura se mantiene constante como consecuencia directa de que la derivada espacial en ese extremo sea nula, es decir, el flujo sea nulo. Para comprobar estas ideas, representamos el flujo en ambos extremos de la varilla:

Representación del flujo en los extremos de la varilla con respecto al tiempo.

Podemos ver al igual que en la otra gráfica como en el extremo derecho no existe flujo y en el extremo izquierdo, el flujo será saliente. Como vimos en la primera gráfica, el flujo empieza en [math] 0 [/math] debido a su condición inicial. Rápidamente aumenta a valores más negativos acorde al aumento de la pendiente. En tiempo [math] t=0.2 ~s[/math] el flujo se reduce rápidamente, hasta que en [math] t=2~s [/math] el flujo desaparece y se llega a la solución estacionaria.

6.1 Simulación

Como hemos hecho en los casos anteriores, modelizamos el cambio de temperatura que se produce en la varilla:

clear all
close all
 
% COEFICIENTES DE FOURIER CALCULADOS MEDIANTE LA APROXIMACIÓN POR LA FÓMULA DEL TRAPECIO
c_k=@(x,k) trapz(x,max(zeros(size(x)),1-4.*abs(x-1/2*ones(size(x)))).*sin((k*pi-pi/2).*x))./trapz(x,(sin((k*pi-pi/2).*x)).^2);
 
% TÉRMINO K-ÉSIMO DE LA SOLUCIÓN DEL PROBLEMA HABIENDO RESTADO LA SOLUCIÓN ESTACIONARIA
w_k=@(x,t,k,D,coef_fourier) coef_fourier.*sin((k*pi-pi/2).*x).*exp(-D.*((-pi/2+k*pi)^2.*t)); 
 
% DEFINIMOS LA SOLUCIÓN ESTACIONARIA
sol_estacionaria=@(x) zeros(size(x));
 
% VALORES DE LOS COEFICIENTES
coeficiente_conductividad_termica=1;
calor_especifico=1;
 
% Definir la longitud y radio de la varilla
longitud = 1;
radio = 0.25;
 
% Definir el número de puntos en la varilla
num_puntos_longitudinal = 1001; % Número de puntos a lo largo de la longitud
num_puntos_circunferencia = 200; % Número de puntos alrededor de la circunferencia
 
% Definir el tiempo y el paso de tiempo
tiempo_final = 2; 
funcion_velocidad_t=@(t) t*1.05;
tiempo=[0,0.001];
while tiempo(end)<tiempo_final
    tiempo(end+1)=funcion_velocidad_t(tiempo(end));
end
 
% Generar la malla de puntos en la varilla
longitudinal = linspace(0, longitud, num_puntos_longitudinal);
circunferencia = linspace(0, 2*pi, num_puntos_circunferencia);
[Longitudinal, Circunferencia] = meshgrid(longitudinal, circunferencia);
 
% Crear un vídeo
pelicula=VideoWriter('Ej1_Ap8_Tuberia','MPEG-4'); %creo el video
pelicula.FrameRate=11; %controla velocidad
open(pelicula); %abrimos el vídeo
 
% Bucle para cada paso de tiempo
for t = tiempo
    figura=figure(1);
    
    % Calcular la temperatura en cada punto de la varilla en el tiempo t
    w = zeros(size(Longitudinal));
    for k = 1:1000
        w = w + w_k(Longitudinal,t,k,coeficiente_conductividad_termica/calor_especifico,c_k(Longitudinal(1,:),k));
    end
    u=w+sol_estacionaria(Longitudinal);
 
    % Graficar la varilla con los colores según la temperatura
    colormap(jet); % Mapa de colores de azul (frío) a rojo (caliente)
    surf(Longitudinal, radio * cos(Circunferencia), radio * sin(Circunferencia), u,'EdgeColor','none');
    colorbar;
    clim([0, 1]); % Fijar el rango de colores
    
    % Añadir etiquetas y ajustes de la figura
    title(['Varilla en t = ' num2str(t)]);
    xlabel('Longitud');
    axis([0 longitud -radio radio -radio radio]);
    view(3);
    
    % Añadir el frame al vídeo
    imagen=getframe(figura);
    writeVideo(pelicula,imagen); %inserto la imagen
    if t==0
        imwrite(imagen.cdata,"Ej1_Ap8_Tuberia_Inicio.png") % guardamos la imagen en el instante inicial
    end
    
end
 
% Cerrar el vídeo
close(pelicula);

Representación la temperatura en la varilla

Esta animación nos afirma todas las conclusiones obtenidas a raíz de las otras gráficas. En concreto, se puede apreciar con claridad como cerca del extremo derecho la temperatura se mantiene constante. Con respecto a los cambios realizados en este código con respecto a las anteriores simulaciones, se han modificado los tiempos para que el video capture más imágenes al principio, donde ocurren cambios significativos con mayor frecuencia, y luego tome imágenes en intervalos más espaciados de tiempo, ya que los cambios son menos notorios.

6.2 Principio del máximo

Por último, vamos a comprobar que el principio del máximo se verifica para este problema, con las condiciones tomadas en el apartado anterior.

Para este último caso, el espacio es [math]Q_{T}=(0,1)\times(0,2) [/math] y se verifica que la solución [math] u(x,t)=\sum_{k=1}^\infty c_{k} sin\left(\left(k \pi-\frac{\pi}{2}\right)x\right)e^{(k \pi -\frac{\pi}{2})^2 Dt }+x :=\sum_{k=1}^\infty u_k (x,t)[/math] converge uniformemente y las derivadas de [math]u_k (x,t)[/math] también. Por tanto, podemos afirmar que:

[math] u_x(x,t) =\sum_{k=1}^\infty \frac{d}{dx} u_k (x,t) [/math]

[math] u_t(x,t) =\sum_{k=1}^\infty \frac{d}{dt}u_k (x,t) [/math]

Como podemos derivar para cualquier orden, se concluye que [math] u(x,t) \in C^{\infty}(Q_{T}) [/math].

Además, [math] u_t(x,t)-u_{xx}(x,t)= \sum_{k=1}^\infty \frac{d}{dt} u_k (x,t)- \frac{d^2}{dx^2}u_k (x,t) =0 \leq 0[/math]. Como las condiciones del principio del máximo se verifican, el máximo de la función ha de estar en la frontera parábolica de [math] Q_T [/math]. En nuestro caso, se puede apreciar como el máximo de la solución es el valor [math]1[/math], que se alcanza en el tiempo inicial[math]t=0 s[/math] y espacio [math]x=0.5 m[/math], perteneciendo claramente a la frontera parábolica de nuestro problema.

7 Solución fundamental de la ecuación del calor

En esta sección vamos a estudiar una solución al sistema

[math] \begin{cases} u_t(x,t)-u_{xx}=0 & x \in \mathbb{R}, t\gt0\\ u(x,0)=f(x) & x \in \mathbb{R} \end{cases} [/math]

Para ello, vamos a usar la solución fundamental de la ecuación del calor. Cabe destacar que como este sistema no tiene condiciones frontera, pues [math] x \in \mathbb{R} [/math], es posible que no haya unicidad de las soluciones.

7.1 Definición solución fundamental

La solución fundamental de la ecuación del calor viene dada por

[math] \Gamma_D (x,t) = \frac{1}{(4 \pi Dt)^{n/2}}e^{-\frac{|x|^2}{4Dt}}, \quad t\gt0 [/math]

Si el lector lo desea, puede consultar la deducción de esta solución en la página 42 del libro Partial differential equations in action from modelling to theory de Sandro Salsa.

Destacamos que, buscando la conservación de la masa total para cualquier tiempo, es lógico que se verifique que [math]\int_{\mathbb{R^{n}}} u(x,t) d\mathbf{x} =q[/math] , con [math] q [/math] constante para todo tiempo. Sin embargo, en la deducción de esta fórmula se ha impuesto la condición [math] q=1 [/math], dando como resultado una familia de funciones Gaussianas parametrizadas en el tiempo.

7.2 Caso unidimensional

En concreto, en el caso unidimensional, la solución fundamental de la ecuación del calor es

[math] \Gamma_D (x,t) = \frac{1}{\sqrt{4 D \pi t}}e^{-\frac{|x|^2}{4Dt}}, \quad t\gt0 [/math].

Tomando el valor de la constante de difusión [math] D=1 [/math], [math] x \in [-1,1] [/math] y [math] t \in [10^{-2},1] [/math]; hemos elaborado el siguiente código de Matlab para representar la solución

clear all
close all

%DEFINIMOS LA SOLUCIÓN FUNDAMENTAL
sol_fundamental=@(x,t,D) 1./(sqrt(4.*pi.*D.*t)).*exp((-x.^2)./(4.*D.*t));

%VALORES DE LOS COEFICIENTES
coeficiente_conductividad_termica=1;
conduccion_calor=1;

%DATOS PARA LA DISCRETIZACIÓN DEL ESPACIO
extremo_izquierdo=-1;
extremo_derecho=1;
paso_espacio=0.001;

%DATOS PARA LA DISCRETIZACIÓN DEL TIEMPO
tiempo_inicial=10^(-2);
tiempo_final=1;
paso_tiempo=0.001;

%DISCRETIZACIÓN DEL ESPACIO Y DEL TIEMPO
[X,T]=meshgrid(extremo_izquierdo:paso_espacio:extremo_derecho,tiempo_inicial:paso_tiempo:tiempo_final); % matrices de mallado espacio-temporal

%REPRESENTAMOS LA SOLUCIÓN FUNDAMENTAL
mesh(X,T,sol_fundamental(X,T,coeficiente_conductividad_termica/conduccion_calor))

Representación de la solución fundamental en una dimensión para [math] x \in [-1,1] [/math] y [math] t \in [10^{-2},1] [/math] .

En este código, lo primero que definimos es la solución fundamental en una dimensión. A continuación, simplemente añadimos los valores del coeficiente de conductividad térmica y del calor específico, que en este caso ambos son iguales a [math] 1 [/math]. Además, al igual que sucedía en los códigos anteriormente explicados, en Matlab trabajamos de manera discreta. Por ello, lo siguiente a definir son los datos que vamos a necesitar para la discretización del espacio y del tiempo y el mallado espaciotemporal haciendo uso de la función de Matlab meshgrid, para así poder hacer la representación gráfica de la solución con respecto al espacio y al tiempo con la función de Matlab mesh.

En la gráfica obtenida como resultado, se aprecian claramente algunas propiedades. La primera de ellas es la paridad de la función en la variable [math] x [/math]. De hecho, de manera analítica se puede comprobar que:

[math] \Gamma_D (x,t) = \frac{1}{(4 \pi Dt)^{n/2}}e^{-\frac{|x|^2}{4Dt}}=\frac{1}{(4 \pi Dt)^{n/2}}e^{-\frac{|-x|^2}{4Dt}}=\Gamma_D (-x,t) [/math]

Demostrándose así la simetría. Además, lo más destacable de la gráfica es el fenómeno que se produce cuando nos aproximamos a [math] t=0 [/math]. Para poder realizar un mejor estudio de este fenómeno, hemos representado la solución fundamental para [math] x \in [-10,10] [/math] y para [math] t \in [10^{-2},5] [/math].

Representación de la solución fundamental en una dimensión para [math] x \in [-10,10] [/math] y [math] t \in [10^{-2},5] [/math] .

En esta gráfica, se aprecia perfectamente que cuando [math] t [/math] tiende a [math] 0 [/math], la solución fundamental tiende a [math] 0 [/math] en cualquier [math] x \neq 0 [/math] y tiende a [math] \infty [/math] en [math] x=0 [/math]. Esta distribución se conoce como la Delta de Dirac.

Se define la Delta de Dirac como el límite de una sucesión de funciones [math] \{f^{n}\}_{n \in \mathbb{N}} [/math] que tienden a cero en todo punto del espacio excepto en un punto para el cual diverge a infinito. Concretamente, llamaremos Delta de Dirac al siguiente límite: [math] \delta = \lim_{n \to \infty} f^n(x) [/math] con [math] \lim_{n \to \infty} f^n(0)=\infty [/math] y [math] \lim _{n \to \infty} f^n(x)=0 [/math] para [math] x \neq 0[/math] que además verifica [math] \int_{\mathbb{R}} f^n(x) \cdot dx =1 [/math].

Además, vemos como el máximo de la solución fundamental se alcanza en [math] (x,t)=(0,0) [/math], cumpliéndose así el principio del máximo. Además, el valor de la función disminuye rápidamente a medida que nos alejamos de este punto en ambas direcciones. Finalmente, cabe destacar que la solución esatacionaria correspondiente a este problema es [math] v(x)=0 [/math]. De hecho,

[math] \lim_{t \to \infty} \Gamma_1= \lim_{t \to \infty} \frac{1}{(4 \pi t)^{n/2}}e^{-\frac{|x|^2}{4t}}=0 [/math]

Sin embargo, la solución fundamental nunca se anula si [math] t\gt0 [/math]. Esto quiere que el sistema tiende a la solución estacionaria, pero no la alcanza en tiempo finito.

7.2.1 Simulación

A continuación, se presenta un código correspondiente a la simulación de una varilla metálica sobre la que se difunde el calor según la solución fundamental. El funcionamiento y construcción de este código es análogo a los ya mencionados en otras secciones. El código sería el siguiente:

clear all
close all
 
% Definimos la solución fundamental
sol_fundamental=@(x,t,D) 1./(sqrt(4.*pi.*D.*t)).*exp((-x.^2)./(4.*D.*t));
 
% Valores de los coeficientes
coeficiente_conductividad_termica=1;
calor_especifico=1;
 
% Definir la longitud y radio de la varilla
longitud_positiva = 1;
radio = 0.25;
 
% Definir el número de puntos en la varilla
num_puntos_longitudinal = 1001; % Número de puntos a lo largo de la longitud
num_puntos_circunferencia = 200; % Número de puntos alrededor de la circunferencia
 
% Definir el tiempo y el paso de tiempo
tiempo_final = 2; 
paso_tiempo = 0.005;
tiempo = 10^(-3):paso_tiempo:tiempo_final;
 
% Generar la malla de puntos en la varilla
longitudinal = linspace(-longitud_positiva, longitud_positiva, num_puntos_longitudinal);
circunferencia = linspace(0, 2*pi, num_puntos_circunferencia);
[Longitudinal, Circunferencia] = meshgrid(longitudinal, circunferencia);
 
% Crear un vídeo
pelicula=VideoWriter('Ej2_Ap1_Tuberia','MPEG-4'); %creo el video
pelicula.FrameRate=27; %controla velocidad
open(pelicula); %abrimos el vídeo
 
% Bucle para cada paso de tiempo
for t = tiempo
    figura=figure(1);
    
    % Calcular la temperatura en cada punto de la varilla en el tiempo t
    u=sol_fundamental(Longitudinal,t,coeficiente_conductividad_termica/calor_especifico);
 
    % Graficar la varilla con los colores según la temperatura
    colormap(jet); % Mapa de colores de azul (frío) a rojo (caliente)
    surf(Longitudinal, radio * cos(Circunferencia), radio * sin(Circunferencia), u,'EdgeColor','none');
    colorbar;
    clim([0, 1]); % Fijar el rango de colores
    
    % Añadir etiquetas y ajustes de la figura
    title(['Varilla con temperatura en t = ' num2str(t)]);
    xlabel('Espacio');
    axis([-longitud_positiva longitud_positiva -radio radio -radio radio]);
    view(3);
    
    % Añadir el frame al vídeo
    imagen=getframe(figura);
    writeVideo(pelicula,imagen); %inserto la imagen
    if t==0
        imwrite(imagen.cdata,"Ej2_Ap1_Tuberia_Inicio.png") % guardamos la imagen en el instante inicial
    end
    
end
 
% Cerrar el vídeo
close(pelicula);

Representación de la varilla metálica con la solución fundamental.

Este caso se da en la vida real cuando el radio de la varilla es despreciable en comparación con la longitud de esta. Sin embargo, en esta simulación, hemos representado una varilla de longitud [math] 2m [/math] para que coincidan los resultados con la primera gráfica de la solución fundamental expuesta.

Como se puede observar, el tiempo inicial del video es [math] t=10^{-2} ~s[/math], puesto que la solución fundamental tiende a la delta de Dirac cuando [math] t \rightarrow 0 [/math], y por lo tanto, la varilla tendría temperatura infinita en sus puntos centrales y cero en el resto de puntos. Por otro lado, se puede ver como las temperaturas más elevadas inicialmente se distribuyen en el centro de la varilla y con el tiempo se van difundiendo hacia los extremos tendiendo la temperatura de la varilla a [math] 0 ºC [/math], aunque nunca alcanzando el estado estacionario.

7.3 Resolución de la ecuación del calor con la solución fundamental

En esta sección, buscamos resolver el siguiente problema:

[math]\left \{ \begin{array}{ll} u_t-u_{xx}=0& \quad x\gt0, t\gt0 \\ u(0,t)=1 & \quad t\gt0 \\ u(x,t)=0 & \quad t\gt0 , x \rightarrow \infty \\ u(x,0)=0 & \quad x\gt0 \end{array} \right. [/math]

Llamamos soluciones autosemejantes a aquellas para las que existe [math] \alpha \in \mathbb{R} [/math] tal que [math] u(\lambda x, \lambda^2 t)= \lambda ^{\alpha} u(x,t) \quad \forall \lambda \gt 0, \forall (x,t) \in Q_{T} [/math].

En primer lugar, vamos a buscar las soluciones autosemejantes de la ecuación del calor [math] u_t-u_{xx}=0 [/math] de la forma [math] u(x,t)=U\left( \frac{x}{\sqrt{t}} \right) [/math]. De esta manera, se tiene que

[math] u_t(x,t)=-\frac{x}{2 t^{3 / 2}} U’\left( \frac{x}{\sqrt{t}} \right) [/math] [math] u_{xx}(x,t)= \frac{1}{t} U’’\left( \frac{x}{\sqrt{t}} \right)[/math]

Y, por lo tanto,

[math] u_t(x,t)-u_{xx}=0 \Rightarrow -\frac{x}{2 t^{3 / 2}} U’\left( \frac{x}{\sqrt{t}} \right) - \frac{1}{t} U’’\left( \frac{x}{\sqrt{t}} \right)=0[/math]

Llamando [math] \xi=\frac{x}{\sqrt{t}} [/math], podemos simplificar al expresión

[math] -\frac{1}{2t} \xi U’\left( \xi \right) - \frac{1}{t} U’’\left( \xi \right) =0 \Rightarrow \frac{1}{t} \left[ \frac{1}{2} \xi U’\left( \xi \right) + U’’\left( \xi \right) \right] =0[/math]

Definiendo [math] f(\xi)= U’\left( \xi \right) [/math] y aplicando el método de separación de variables, se llega a que

[math] \frac{df}{d \xi}=-\frac{1}{2} \xi f \Rightarrow \int \frac{df}{f}=-\frac{1}{2} \int \xi d\xi \Rightarrow U’\left( \xi \right) = f(\xi)=D e^{-\xi^2 /4}[/math] , con [math] D \in \mathbb{R} [/math]

A continuación, aplicando el teorema fundamental del cálculo, se deduce que

[math] U(\xi)=D \int_{0}^{\xi} e^{-\frac{\eta^2}{4}} d\eta +E, \quad D,E \in \mathbb{R} [/math]

Para calcular los valores de [math] D,E \in \mathbb{R} [/math], debemos hacer uso de las condiciones iniciales y de frontera de las que disponemos. En primer lugar, se debe cumplir que [math] u(0,t)=1 [/math] o, lo que es lo mismo,

[math] U(0)= D \int_{0}^{0} e^{-\frac{\eta^2}{4}} d\eta +E=E=1 [/math]

Por lo tanto, [math] U(\xi)=D \int_{0}^{\xi} e^{-\frac{\eta^2}{4}} d\eta +1, \quad D \in \mathbb{R} [/math]. Por otra parte, el ejercicio exige que [math] u(x,t) \rightarrow 0 [/math] cuando [math] x \rightarrow \infty [/math]. Esto es equivalente a que [math] U(\xi) \rightarrow 0 [/math] cuando [math]\xi \rightarrow \infty [/math] y haciendo uso de la paridad de la función [math] e^{-\frac{\eta^2}{4}}, ~\eta \in \mathbb{R} [/math], se obtiene la siguiente serie de igualdades:

[math]0= \lim_{\xi \to \infty } U(\xi)= \lim_{\xi \to \infty} D \int_{0}^{\xi} e^{-\frac{\eta^2}{4}} d\eta + 1= \frac{D}{2} \int_{\mathbb{R}} e^{-\frac{\eta^2}{4}} d\eta +1 [/math]

Además, si aplicamos el cambio de variable [math] z=\frac{\eta}{2} [/math] y sabiendo que [math] \int_{\mathbb{R}} e^{-z^2} dz= \sqrt{\pi} [/math], llegamos a que

[math] D \sqrt{\pi}+1=0 \Rightarrow D=-\frac{1}{\sqrt{\pi}} [/math]

De esta manera, concluimos que la solución buscada es

[math] U(\xi)=-\frac{1}{\sqrt{\pi}} \int_{0}^{\xi} e^{-\frac{\eta^2}{4}} d\eta +1=1-\frac{2}{\sqrt{\pi}} \int_{0}^{\xi/2} e^{-z^2} dz=1-erf\left(\frac{\xi}{2}\right) [/math]

siendo [math] erf(x)= \frac{2}{\sqrt{\pi}} \int_{0}^{x} e^{-z^2} dz [/math].

Por lo tanto, deshaciendo el cambio de variable,

[math] u(x,t)=U\left(\frac{x}{\sqrt{t}}\right)=1-erf\left(\frac{x}{2\sqrt{t}}\right) [/math]

Cabe destacar que esta solución también verifica la condición inicial que faltaba por verificar, pues [math] \lim_{t \to 0} u(x,t)= \lim_{t \to 0} 1- \frac{2}{\sqrt{\pi}}\int_{0}^{\frac{x}{2 \sqrt{t}}} e^{-z^2} dz=1- \frac{2}{\sqrt{\pi}}\frac{1}{2}\int_{\mathbb{R}} e^{-z^2} dz=1- \frac{2}{\sqrt{\pi}}\frac{1}{2} \sqrt{\pi}=0 [/math].

Además,

[math] \lim_{t \to \infty} u(x,t)= \lim_{t \to \infty} 1-erf\left(\frac{x}{2\sqrt{t}}\right)= 1 [/math]

Sin embargo, de nuevo ocurre que la función solución nunca toma el valor 1, pues [math] erf\left(\frac{x}{2\sqrt{t}}\right) [/math] nunca se anula. Esto indica que a pesar de que el sistema tienda a la solución estacionaria <v(x)=1>, esta no se alcanza en tiempo finito.

Una vez calculada la solución al problema, hemos elaborado un código en Matlab para representarla:

clear all
close all
% SOLUCIÓN DEL EJERCICIO 10b
sol_erfc=@(x,t) 1-erf(x./(2.*sqrt(t)));
% DATOS PARA LA DISCRETIZACIÓN DEL ESPACIO
extremo_izquierdo=0;
extremo_derecho=20;
paso_espacio=0.1;
% DATOS PARA LA DISCRETIZACIÓN DEL TIEMPO
tiempo_inicial=0;
tiempo_final=10;
paso_tiempo=0.001;
% DISCRETIZACIÓN DEL ESPACIO Y DEL TIEMPO
[X,T]=meshgrid(extremo_izquierdo:paso_espacio:extremo_derecho,tiempo_inicial:paso_tiempo:tiempo_final); % matrices de mallado espacio-temporal
% REPRESENTAMOS LA SOLUCIÓN
figure(1)
mesh(X,T,sol_erfc(X,T))
title('Representación de la solución')
xlabel('Espacio')
ylabel('Tiempo')
zlabel('Temperatura')

Representación de la solución obtenida de la forma [math] u(x,t)=U\left( \frac{x}{\sqrt{t}} \right) [/math] con [math]x \in (0,20)[/math] y [math] t \in (0,10) [/math].

En este código, primero hemos definido la solución del ejercicio como una función. A continuación, como en el resto de los programas, hemos definido los datos necesarios para poder hacer una discretización del tiempo y del espacio. En este caso, se ha representado la función en [math] (x,t) \in [0,20] \times [0,10] [/math]. Con respecto a la gráfica obtenida como resultado, se puede apreciar a la perfección que se cumple que [math] u(0,t)=1, t\gt0 [/math] y [math] u(x,t) \rightarrow 0 [/math] cuando [math] x \rightarrow \infty [/math]. Además, [math] u(x,0)=0, x\gt0 [/math]. En esta gráfica, es importante destacar que se produce una discontinuidad en [math] (x,t)=(0,0)[/math]. Esta discontinuidad se debe a que las propias condiciones iniciales y de frontera no coinciden en dicho punto. A pesar de esto, Matlab hace una representación continua de la función.

7.3.1 Simulación

De nuevo, hemos hecho la simulación en una varilla metálica para poder ver visualmente cómo se comporta el calor. Como en este caso la solución está definida para todo [math] \forall x\gt0 [/math], hemos representado un trozo de la varilla metálica con una longitud de [math] 100m [/math] con un radio [math] 0.25m [/math].

clear all
close all
 
% SOLUCIÓN DEL EJERCICIO 10b
sol=@(x,t) 1-erf(x./(2.*sqrt(t)));
 
% Valores de los coeficientes
coeficiente_conductividad_termica=1;
calor_especifico=1;
 
% Definir la longitud y radio de la varilla
longitud = 100;
radio = 0.25;
 
% Definir el número de puntos en la varilla
num_puntos_longitudinal = 5001; % Número de puntos a lo largo de la longitud
num_puntos_circunferencia = 200; % Número de puntos alrededor de la circunferencia
 
% Definir el tiempo y el paso de tiempo
tiempo_final = 500; 
paso_tiempo = 1;
tiempo = 0:paso_tiempo:tiempo_final;
 
% Generar la malla de puntos en la varilla
longitudinal = linspace(0, longitud, num_puntos_longitudinal);
circunferencia = linspace(0, 2*pi, num_puntos_circunferencia);
[Longitudinal, Circunferencia] = meshgrid(longitudinal, circunferencia);
 
% Crear un vídeo
pelicula=VideoWriter('Ej2_Ap2_Tuberia','MPEG-4'); %creo el video
pelicula.FrameRate=34; %controla velocidad
open(pelicula); %abrimos el vídeo
 
% Bucle para cada paso de tiempo
for t = tiempo
    figura=figure(1);
    
    % Calcular la temperatura en cada punto de la varilla en el tiempo t
    u=sol(Longitudinal,t);
 
    % Graficar la varilla con los colores según la temperatura
    colormap(jet); % Mapa de colores de azul (frío) a rojo (caliente)
    surf(Longitudinal, radio * cos(Circunferencia), radio * sin(Circunferencia), u,'EdgeColor','none');
    colorbar;
    clim([0, 1]); % Fijar el rango de colores
    
    % Añadir etiquetas y ajustes de la figura
    title(['Varilla con temperatura en t = ' num2str(t)]);
    xlabel('Espacio');
    axis([0 longitud -radio radio -radio radio]);
    view(3);
    
    % Añadir el frame al vídeo
    imagen=getframe(figura);
    writeVideo(pelicula,imagen); %inserto la imagen
    if t==0
        imwrite(imagen.cdata,"Ej2_Ap2_Tuberia_Inicio.png") % guardamos la imagen en el instante inicial
    end
    
end
 
% Cerrar el vídeo
close(pelicula);

Evolución de la solución en una varilla de longitud [math] 100~m [/math] para tiempos [math] t [/math] comprendidos entre [math] 0 ~s [/math] y [math] 500 ~s[/math]

Este caso es muy distinto a los estudiados con anterioridad. El calor se concentra en el extremo izquierdo de la varilla en el instante inicial y con el paso del tiempo, se desplaza hacia la derecha constantemente. Recordamos que el radio de la varilla es despreciable en comparación con su longitud. Si tenemos en cuenta esto, podemos considerar dicha longitud infinita y, por ello, la solución estacionaria no se alcanza en tiempo finito. De hecho, destacamos que el vídeo recién introducido alcanza el tiempo [math] t=500~s[/math] y en dicho instante el calor se ha difundido muy poco en comparación con la longitud de la varilla.

7.4 Solución de la ecuación del calor empleando la convolución

Planteamos ahora un nuevo escenario. En este caso, partimos de la ecuación del calor [math] u_t-u_{xx}=0 [/math] en [math] x \in \mathbb{R} [/math] asociada al dato inicial [math] u_0(x)=1_{[-1,1]} [/math].

En este caso, la solución de dicha ecuación viene dada por la convolución:

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

Como la condición inicial es [math] u_0(x)=1_{[-1,1]} [/math], la integral de la convolución de la solución fundamental y la condición inicial será nula fuera del dominio de integración [math] [-1,1] [/math], por lo que es más sencillo reducir el cálculo de dicha convolución en todo [math] \mathbb{R} [/math] a únicamente el intervalo [math] [-1,1] [/math].

De esta manera, podemos dibujar la solución [math] u(x,t) [/math] en diferentes instantes de tiempo. Para ello, hemos implementado el siguiente código en Matlab:

clear all
close all
 
% DEFINIMOS LA SOLUCIÓN FUNDAMENTAL
sol_fundamental=@(x,t,D) 1./(sqrt(4.*pi.*D.*t)).*exp((-x.^2)./(4.*D.*t));
 
% DEFINIMOS LA FUNCIÓN SOLUCIÓN EN EL INSTANTE INICIAL
u_0=@(x) 1*((-1<=x)&(x<=1));
 
% VALORES DE LOS COEFICIENTES
coeficiente_conductividad_termica=1;
calor_especifico=1;
 
% DEFINIMOS Y COMO LA VARIABLE SOBRE LA QUE VAMOS A INTEGRAR
Y=linspace(-1,1,100);
 
% DATOS PARA LA DISCRETIZACIÓN DEL ESPACIO
extremo_izquierdo=-5;
extremo_derecho=5;
paso_espacio=0.001;
X=[extremo_izquierdo:paso_espacio:extremo_derecho];
 
% TIEMPOS A ESTUDIAR
T=[0.001,0.01,0.1];
t_final=100000;
 
% DEFINIMOS LA SOLUCIÓN u(x,t)
figure(1)
u=zeros(length(T),length(X));
for j=1:length(T)
    for i=1:length(X)
        u(j,i)=trapz(Y,sol_fundamental(X(i)-Y,T(j),coeficiente_conductividad_termica/calor_especifico).*u_0(Y));
    end
    subplot(1,length(T),j)
    plot(X,u(j,:),'b','LineWidth',1.5)
    axis([-2,2,0,1.001])
    xlabel('Espacio')
    ylabel('Temperatura')
    title("t="+num2str(T(j)))
end
sgtitle('Representación de la solución para distintos instantes de tiempo t')
 
% Crear un vídeo
pelicula=VideoWriter('Solucion_reales','MPEG-4'); %creo el video
pelicula.FrameRate=15; %controla velocidad
open(pelicula); %lo abro
 
funcion_velocidad_t=@(t) t*1.1;
vector_t=[0,0.001];
while vector_t(end)<t_final
    vector_t(end+1)=funcion_velocidad_t(vector_t(end));
end
 
for j=1:length(vector_t)
    figura=figure(2);
    for i=1:length(X)
        u(j,i)=trapz(Y,sol_fundamental(X(i)-Y,vector_t(j),coeficiente_conductividad_termica/calor_especifico).*u_0(Y));
    end
    plot(X,u(j,:),'b','LineWidth',1.5)
    axis([extremo_izquierdo,extremo_derecho,0,1.001])
    view(2);
    xlabel('Espacio')
    ylabel('Temperatura')
    title("t="+num2str(vector_t(j)))
    
    % Añadir el frame al vídeo
    imagen=getframe(figura);
    writeVideo(pelicula,imagen); %inserto la imagen
end
 
close(pelicula);

Primero de todo, se define la solución fundamental y la condición inicial, es decir, la función [math] u_0(x) [/math] y los valores de los coeficientes de conductividad térmica y de calor específico.

Se define la variable [math] Y [/math] como la variable sobre la que vamos a integrar. A continuación, se discretiza el espacio definiendo por tanto los extremos donde está definida la solución y el paso de discretización deseado entre dichos puntos. Por último, antes de definir la solución se establece un vector de tiempos [math] T [/math] donde se quiere realizar un estudio de la solución obtenida en dichos instantes de tiempo.

Seguidamente, se define la solución dada por la convolución de las dos funciones ya mencionadas y se aproxima haciendo uso de la fórmula del trapecio en el intervalo espacial definido. En último lugar, se definen los comandos adecuados requeridos para la realización de un video que muestre las soluciones de la ecuación en distintos instantes de tiempo. Concretamente, desde [math]t= 0.0011 ~s [/math] a [math] t=97451 ~s[/math].

En relación con los resultados obtenidos al implementar el código, son destacables varios aspectos.

Representación de la solución [math] u(x,t) [/math] dada por la convolución de la solución fundamental y la condición inicial [math] u_0(x) [/math] para los instantes de tiempo [math] t=0.001~s, ~ t=0.01~s [/math] y [math] t=0.1~s[/math].

En primer lugar, como era de esperar, para valores de [math] t [/math] muy cercanos al valor [math] t=0 [/math] la solución [math] u(x,t) [/math] obtenida es muy similar a la condición inicial [math] u_0(x) [/math].

Evolución de la solución para tiempos [math] t [/math] comprendidos entre [math] t=0 ~s [/math] y [math] t=97451.4 ~s[/math]

En segundo lugar, se puede apreciar como a medida que aumenta el tiempo, la solución [math] u(x,t) [/math] tiende a cero, la solución estacionaria de la ecuación, aunque no la alcanza en tiempo finito. Conceptualmente esto quiere decir que para tiempos lo suficientemente grandes se produce una variación despreciable de la temperatura en el objeto de estudio.

7.4.1 Simulación

Calculada la nueva solución, se presenta el código correspondiente a la simulación de la difusión del calor sobre la varilla. Su construcción es análoga a los anteriores.

clear all
close all
 
% DEFINIMOS LA SOLUCIÓN FUNDAMENTAL
sol_fundamental=@(x,t,D) 1./(sqrt(4.*pi.*D.*t)).*exp((-x.^2)./(4.*D.*t));
 
% DEFINIMOS LA FUNCIÓN SOLUCIÓN EN EL INSTANTE INICIAL
u_0=@(x) 1*((-1<=x)&(x<=1));
 
% Valores de los coeficientes
coeficiente_conductividad_termica=1;
calor_especifico=1;
 
% Definir la longitud y radio de la varilla
longitud = 1;
radio = 0.25;
 
% DEFINIMOS Y COMO LA VARIABLE SOBRE LA QUE VAMOS A INTEGRAR
Y=linspace(-1,1,100);
 
% DATOS PARA LA DISCRETIZACIÓN DEL ESPACIO
extremo_izquierdo=-longitud;
extremo_derecho=longitud;
paso_espacio=0.001;
X=[extremo_izquierdo:paso_espacio:extremo_derecho];
 
% Definir el tiempo y el paso de tiempo
tiempo_final = 10; 
funcion_velocidad_t=@(t) t*1.05;
tiempo=[0.001];
while tiempo(end)<tiempo_final
    tiempo(end+1)=funcion_velocidad_t(tiempo(end));
end
 
% DEFINIMOS Y COMO LA VARIABLE SOBRE LA QUE VAMOS A INTEGRAR
Y=linspace(-1,1,100);
 
% Generar la malla de puntos en la varilla
num_puntos_circunferencia = 200; % Número de puntos alrededor de la circunferencia
circunferencia = linspace(0, 2*pi, num_puntos_circunferencia);
[Longitudinal, Circunferencia] = meshgrid(X, circunferencia);
 
% Crear un vídeo
pelicula=VideoWriter('Ej2_Ap3_Tuberia','MPEG-4'); %creo el video
pelicula.FrameRate=13; %controla velocidad
open(pelicula); %abrimos el vídeo
 
% Bucle para cada paso de tiempo
for t = tiempo
    figura=figure(1);
 
    u=zeros(size(Longitudinal));
 
    % Calcular la temperatura en cada punto de la varilla en el tiempo t
    for i=1:length(Longitudinal(1,:))
        u(1,i)=trapz(Y,sol_fundamental(X(i)-Y,t,coeficiente_conductividad_termica/calor_especifico).*u_0(Y));
    end
    for j=1:length(Longitudinal(:,1))
        u(j,:)=u(1,:);
    end
    % Graficar la varilla con los colores según la temperatura
    colormap(jet); % Mapa de colores de azul (frío) a rojo (caliente)
    surf(Longitudinal, radio * cos(Circunferencia), radio * sin(Circunferencia), u,'EdgeColor','none');
    colorbar;
    clim([0, 1]); % Fijar el rango de colores
    
    % Añadir etiquetas y ajustes de la figura
    title(['Varilla con temperatura en t = ' num2str(t)]);
    xlabel('Espacio');
    axis([-longitud longitud -radio radio -radio radio]);
    view(3);
    
    % Añadir el frame al vídeo
    imagen=getframe(figura);
    writeVideo(pelicula,imagen); %inserto la imagen
    if t==0
        imwrite(imagen.cdata,"Ej2_Ap3_Tuberia_Inicio.png") % guardamos la imagen en el instante inicial
    end
    
end
 
% Cerrar el vídeo
close(pelicula);

Como se puede observar en el vídeo, la varilla metálica está definida en el dominio espacial [math] [-1,1] [/math]. En el instante [math] t=0 [/math] se verifica la condición inicial y, por ello, toda la varilla comienza con una temperatura de [math] 1ºC [/math]. Con el transcurso del tiempo, el calor comienza a difundirse por el resto de la varilla. Es destacable que se ha establecido el tiempo de ejecución de manera que se reproduce mucho más lento y detallado los instantes de tiempo iniciales y, de una manera más rápida a medida que pasa el tiempo. Por último, la temperatura de la varilla tiende a estabilizarse a una temperatura nula.

7.5 Solución fundamental de la ecuación del calor en dimensión 2

Una vez hemos estudiado y analizado cómo se comporta la solución fundamental del calor en una dimensión vamos a analizar brevemente el caso de la solución fundamental en dimensión [math] 2 [/math]. La solución, tal y como hemos mencionado previamente, sería:

[math] \Gamma_D (x,t) = \frac{1}{(4 \pi Dt)^{n/2}}e^{-\frac{|x|^2}{4Dt}} ~t\gt0 [/math]

En dos dimensiones, la geometría del problema puede ser más representativa, especialmente en sistemas bidimensionales como placas o láminas. Por otro lado, el sistema y soluciones pueden resultar más complejos. Este aspecto gana importancia cuando tenemos casos con propiedades significativas del medio o condiciones de contorno o iniciales más elaboradas.

Asimismo, en dimensión [math] 2 [/math] es posible modelar fenómenos que no se pueden capturar en una dimensión, como la difusión del calor a lo largo de un plano sin necesidad de que el flujo sea uniforme en todos los puntos de una misma sección y que la difusión se produzca únicamente en una dirección fijada.

Por último, cabe destacar que la solución estacionaria de la solución fundamental es:

[math] \lim_{t \to \infty} \Gamma_D (x,t) = \lim_{t \to \infty} \frac{1}{(4 \pi t)^{n/2}}e^{-\frac{|x|^2}{4t}}=0 [/math]

Sin embargo, de nuevo vuelve a ocurrir que la solución fundamental no se anula en ningún punto para t>0. Esto quiere decir que la solución fundamental tiende a la solución estacionaria, pero no la alcanza en tiempo finito.

La manera más sencilla de entenderlo es con una representación gráfica. Por tanto, hemos realizado un código en Matlab que implementa dicha solución en distintos instantes de tiempo. El código es el siguiente:

clear all
close all
 
% DEFINIMOS LA SOLUCIÓN FUNDAMENTAL
sol_fundamental=@(x1,x2,t,D) 1./((4.*pi.*D.*t)).*exp((-sqrt(x1.^2+x2.^2).^2)./(4.*D.*t));
 
% VALORES DE LOS COEFICIENTES
coeficiente_conductividad_termica=1;
calor_especifico=1;
 
% DATOS PARA LA DISCRETIZACIÓN DEL ESPACIO
extremo_izquierdo=-1;
extremo_derecho=1;
paso_espacio=0.001;
[X1,X2]=meshgrid(extremo_izquierdo:paso_espacio:extremo_derecho,extremo_izquierdo:paso_espacio:extremo_derecho);
 
% TIEMPOS A ESTUDIAR
T=[0.00001,0.001,0.01,0.1];
t_final=0.001;
 
% DEFINIMOS LA SOLUCIÓN u(x,t)
figure(1)
for j=1:length(T)
    subplot(1,length(T),j)
    mesh(X1,X2,sol_fundamental(X1,X2,T(j),coeficiente_conductividad_termica/calor_especifico))
 
    xlabel('Valores de x_1')
    ylabel('Valores de x_2')
    zlabel('Temperatura')
    title("t="+num2str(T(j)))
end
sgtitle('Representación de la solución fundamental bidimensional para distintos instantes de tiempo t')
 
% Crear un vídeo
pelicula=VideoWriter('Solucion_fundamental_bidimensional','MPEG-4'); %creo el video
pelicula.FrameRate=10; %controla velocidad
open(pelicula); %lo abro
 
 
funcion_velocidad_t=@(t) t*0.95;
vector_t=[0.05];
while vector_t(end)>t_final
    vector_t(end+1)=funcion_velocidad_t(vector_t(end));
end
 
for j=1:length(vector_t)
    figura=figure(2);
    mesh(X1,X2,sol_fundamental(X1,X2,vector_t(j),coeficiente_conductividad_termica/calor_especifico))
    axis([-1,1,-1,1,0,90])
    view(3);
    xlabel('Valores de x_1')
    ylabel('Valores de x_2')
    zlabel('Temperatura')
    colormap(jet);
    colorbar
    clim([0,90])
 
    title("t="+num2str(vector_t(j)))
 
    % Añadir el frame al vídeo
    imagen=getframe(figura);
    writeVideo(pelicula,imagen); %inserto la imagen
end
 
close(pelicula);

La construcción de este código es análoga a la ya realizada en la sección 7.2. La única diferencia sería la definición de la solución fundamental, ya que habría que introducir la correspondiente a dimensión [math] 2 [/math], además de ajustar la definición de los mallados, que ahora sería un mallada únicamente espacial. Además, se hace una gráfica en 3 dimensiones para cada valor del tiempo.

Representación de la solución fundamental de la ecuación del calor [math] \Gamma_D (x,t) [/math] en dimensión [math] n=2 [/math] para los instantes de tiempo [math] t=0.00001~s, ~ t=0.001~s,~t=0.01~s [/math] y [math] t=0.1[/math].

Como conclusión de estas representaciones se puede observar que cuanto más cerca es el valor de la [math] t [/math] a cero, más singular es la solución. Este fenómeno se aprecia a la perfección en el siguiente vídeo, en el que el tiempo se va aproximando a [math] 0 ~s [/math]:

Representación de la solución fundamental de la ecuación del calor bidimensional con una sucesión decreciente de tiempos

Concretamente, al igual que en el caso unidimensional, la solución fundamental de la ecuación del calor en dimensión [math] n=2 [/math] coincide con la delta de Dirac (bidimensional) para [math] t=0 [/math]. Esta observación se aprecia perfectamente en el vídeo realizado. La delta de Dirac en dos dimensiones se define de manera similar a la delta de Dirac unidimensional:

[math] \delta(x, y) = \begin{cases} \infty & \text{si } (x, y) = (0, 0) \\ 0 & \text{en cualquier otro lugar del plano} \end{cases}[/math]

Por otra parte, recordando que en el vídeo el tiempo decrece, se aprecia claramente cómo tiende a la solución estacionaria, [math] v(\mathbf{x})=0 [/math] cuando [math] t [/math] crece.

8 Referencias

SectionB. Similarity solutions.
Ecuaciones en derivadas parciales Facultad de Matemáticas. Universidad de Sevilla
Mauro Milonga Miguel. Trabajo de fin de grado: Ecuación del calor. Universidad de Murcia
[Partial differential equations in action from modelling to theory. Sandro Salsa]

@@ Línea 591: / Línea 591: @@
 La construcción del código es análoga a la ya mencionada en la sección 3.3, pero con una pequeña diferencia. Se define un vector de coeficientes de conductividad térmica <math> k </math> y se repite el proceso de representación para cada uno de los valores.
-[[Archivo: Ej1 Ap6 Tuberia2.gif|600px|thumb|center| Simulación de la solución en la varilla metálica. Nota: Si no se visualiza la animación de este archivo gif, pinchar en el archivo para su visualización.]]
+[[Archivo: Ej1 Ap6 Tuberia2.gif|900px|thumb|center| Simulación de la solución en la varilla metálica. Nota: Si no se visualiza la animación de este archivo gif, pinchar en el archivo para su visualización.]]
 En este caso, se aprecian perfectamente las observaciones realizadas. El coeficiente de conductividad térmica afecta principalmente a la velocidad con la que se transmite el calor. Si observamos el primer y tercer video, la velocidad a la que se transmite el calor es muy distinta. En el primer caso, con coeficiente de difusión <math> D=1 </math>, el sistema alcanza el estado estacionario mucho más rápido. De hecho, el caso en el que el coeficiente de difusión <math> D=1/6 </math>, el sistema aún no ha alcanzado el estado estacionario cuando finaliza el vídeo.
@@ Línea 1231: / Línea 1231: @@
 }}
-[[Archivo:Ej2 Ap1 Tuberia.gif|600px|thumb|center| Representación de la varilla metálica con la solución fundamental.]]
+[[Archivo: Ej2 Ap1 Tuberia.gif|600px|thumb|center| Representación de la varilla metálica con la solución fundamental.]]
 Este caso se da en la vida real cuando el radio de la varilla es despreciable en comparación con la longitud de esta. Sin embargo, en esta simulación, hemos representado una varilla de longitud <math> 2m </math> para que coincidan los resultados con la primera gráfica de la solución fundamental expuesta.
@@ Línea 1693: / Línea 1693: @@
 * [https://webs.um.es/gustavo.garrigos/tfg/Milonga_TFG_sept2020.pdf. Mauro Milonga Miguel. Trabajo de fin de grado: Ecuación del calor. Universidad de Murcia]
 *[Partial differential equations in action from modelling to theory. Sandro Salsa]
+[[Categoría:EDP]]
+[[Categoría:EDP23/24]]