Ecuación del calor (Raúl, Sofía, Jaime)
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Ecuación del calor |
| Asignatura | EDP |
| Curso | 2023-24 |
| Autores | Raúl Ortega
Sofía Gómez Jaime Sáenz de Miera |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
En este artículo vamos a estudiar la ecuación del calor en una dimensión mediante el estudio de una varilla metálica de longitud 1 y que se encuentra aislada por su superficie lateral, es decir, la conducción de calor solo se produce en la dirección longitudinal. Además, vamos a ver cómo varía la solución en función de la variación de distintos parámatros como las condiciones iniciales o el coeficiente de difusión.
Contenido
1 Preliminares
Antes de empezar a desarrollar el problema necesitamos conocer algunos conceptos que nos van a permitir modelizarlo y sacar sus ecuaciones.
- La ley de Fourier establece que el flujo de transferencia de calor es proporcional y de sentido contrario al gradiente de temperatura en esa dirección.
Siendo k la conductividad térmica.
- Conservación de la energía.
2 Problema original
En este primer apartado, vamos a suponer que la temperatura inicial de la varilla es 0ºC y que en el extremo izquierdo se consigue mantener la temperatura a 0ºC mientras que en el derecho la temperatura es siempre de 1ºC. Teniendo esto en cuenta y los conceptos desarrollados en el apartado anterior, el sistema de ecuaciones que representa el problema es:
[math] \begin{cases} u_t-u_{xx}=0 & x\in(0,1),t\gt0\\ u(0,t)=0&t\gt0\\ u(1,t)=1&t\gt0\\ u(x,0)=0&x\in(0,1) \end{cases} [/math]
2.1 Homogeneización
Antes de empezar con la resolución del problema vamos a homogeneizar las condiciones de contorno. Para hacer esto vamos a restarle al problema la solución estacionaria, es decir, la que permanece invariante respecto al tiempo. Esto implica que [math] u_t\sim 0 [/math] y que [math] u(x,t)\sim v(x) [/math]. Por lo que el problema queda de la siguiente forma:
[math] \begin{cases} v_{xx}=0 & x\in(0,1)\\ v(0)=0\\ v(1)=1\\ \end{cases} [/math]
La solución de este problema, y por tanto la solución estacionaria es: [math] v(x)=x [/math].
v=@(x,t) x;
colormap jet
fsurf(v, [0 1 0 1], 'FaceColor','interp')
xlabel('x')
ylabel('t')
zlabel('v(x)')
Una vez tenemos la solución estacionaria, tomamos [math] w(x,t)=u(x,t)-v(x) [/math] y vemos como queda el problema homogeneizado para [math] w(x,t)[/math]:
[math] \begin{cases} w_t-w_{xx}=0 & x\in(0,1),t\gt0\\ w(0,t)=0&t\gt0\\ w(1,t)=0&t\gt0\\ w(x,0)=-x&x\in(0,1) \end{cases} [/math]
2.2 Resolución
Una vez homogeneizado el sistema podemos empezar a buscar la solución por el método de separación de variables, en el cual vamos a suponer que la solución es de la forma: [math] w(x,t)=X(x)T(t) [/math]. Como [math]w(x,t)[/math] es solución obtenemos que se tiene que cumplir [math]\frac{X''(x)}{X(x)}=\frac{T'(t)}{T(t)}=\lambda_k[/math] dando lugar a la resolución de dos problemas por separado, uno que depende de [math] x [/math] y otro que depende de [math] t[/math]. Al solucionar ambos problemas obtenemos el siguiente resultado: [math] w(x,t)=\sum_{k=1}^\infty c_k\sin(k\pi x)e^{-k^2\pi^2 t} [/math].
Para que [math] w(x,t) [/math] sea solución de nuestro problema, solo falta ver que cumpla la condición incial, es decir, que [math] w(x,0)=c_k\sin(k\pi x)=-x [/math]. De esta igualdad sacamos que si extendemos la función [math] -x [/math] de forma impar los coeficientes [math] c_k[/math] van a ser los coeficiente de la serie de Fourier. Por tanto, la solución queda de la siguiente manera: [math] w(x,t)=\sum_{k=1}^\infty\frac{2(-1)^k}{k\pi}\sin(k\pi x)e^{-k^2\pi^2 t}[/math]
Esta función es la solución del problema homogeneizado, por lo que para obtener la solución del problema original tenemos que deshacer el cambio de variable. La solución es:
Ahora vamos a graficar los 10 primeros términos de la serie en [math] t \in [0,1][/math].
clear all
%Definimos las funciones
wi=@(x,t,k) ((2*(-1)^(k))/(k*pi)).*sin(k.*pi.*x).*exp(-k^2.*pi^2.*t);
v=@(x) x;
%Sumamos los 10 primeros términos
w=@(x,t) 0;
for n=1:10
w=@(x,t) w(x,t) + wi(x,t,n);
end
u=@(x,t) w(x,t)+v(x);
% Gráfica en 3D
close all
colormap jet
fsurf(u,[0 1 0 1], 'FaceColor','interp')
xlabel('x')
ylabel('t')
zlabel('u(x,t)')
En la gráfica podemos observar que la solución del poblema alcanza la solución estacionaria.
2.3 Flujo en los extremos
Por la ley de Fourier sabemos que [math] q=-k u_x [/math], siendo q el flujo. Por lo que para interpretar como es el flujo en nuestro problema, vamos a graficar como es en los extremos suponiendo que [math] k=1 [/math].
%Definimos la función
dwi=@(x,t,k) (2*(-1)^(k+1)).*cos(k.*pi.*x).*exp(-k^2.*pi^2.*t);
dv=@(x) 1;
%Sumamos los 10 primeros términos
dw0=@(t) 0;
dw1=@(t) 0;
for n=1:10
dw0=@(t) dw0(t) + dwi(0,t,n);
dw1=@(t) dw1(t) + dwi(1,t,n);
end
du0=@(t) dw0(t)-dv(0);
du1=@(t) dw1(t)-dv(1);
% Gráfica
close all
t=linspace(0,1,100);
subplot(2,1,1)
plot(t,du0(t),color='b',LineWidth=1.5)
xlabel('t')
ylabel('u´(0,t)')
title('Flujo en el extremo izquierdo')
subplot(2,1,2)
plot(t,du1(t),color='b',LineWidth=1.5)
xlabel('t')
ylabel('u´(1,t)')
title('Flujo en el extremo derecho')
Al observar la gráfica en el extremo izquierdo vemos que al principio no hay flujo debido a que tanto la barra como el extremo están a 0º. Una vez avanza el tiempo el flujo va siendo cada vez más negativo lo que implica que cada vez sale más flujo por el extremo izquierdo. Esto tiene sentido ya que como la barra se va calentando, como podemos ver en la gráfica de la solución, y el extremo se mantiene a 0º el flujo tiende ir de donde hay mayor temperatura a donde hay una temperatura menor.
Por el contrario, en la gráfica del extremos derecho podemos observar que inicialmente hay mucho flujo hacia el interior de la barra debido a que la barra está a 0º y el extremo a 1º. Una vez avanza el tiempo el flujo va siendo cada vez menos negativo lo que implica que cada vez entra menos flujo por el extremo derecho. Esto es coherente ya que como la barra se va calentando, y el extremo se mantiene a 1º el flujo tiende ir de donde hay mayor temperatura a donde hay una temperatura menor como hemos mencionado antes.
Al comparar ambas gráficas, cabe destacar dos aspectos característicos:
- Lo primero que cabe preguntarse es por qué cuando en la gráfica hay una gran variación del flujo, en el extremo izquierdo desciende más moderadamente, y en el extremo derecho asciende drásticamente. Esto se debe a que en el extremo izquierdo la diferencia de temperatura entre dos puntos cercanos es menor que en el derecho, cosa que podemos observar en la gráfica de la solución.
- También llama la atención que ambas gráficas se acaban estabilizando en un flujo de -1, este punto coincide con el punto en el que la solución alcanza la solución estacionaria, por lo que tiene sentido que el flujo en ambos extremos se estabilice en el mismo valor.
3 Variación coeficiente de conductividad
En este apartado vamos a ver como varia la solución cuando variamos el coeficiente de conductividad térmica de 1 a 1/2. Con este coeficiente el sistema queda de la siguiente forma:
[math] \begin{cases} u_{1t}-\frac{1}{2}u_{1xx}=0 & x\in(0,1),t\gt0\\ u_1(0,t)=0&t\gt0\\ u_1(1,t)=1&t\gt0\\ u_1(x,0)=0&x\in(0,1) \end{cases} [/math]
Resolviendolo de la misma forma que el anterior obtenemos que la solución es:
y que la estacionaria no varía [math] v_1(x)=x [/math].
clear all
%Definimos las funciones
wi=@(x,t,k) ((2*(-1)^(k))/(k*pi)).*sin(k.*pi.*x).*exp(-k^2.*pi^2.*t./2);
v=@(x) x;
%Sumamos los 10 primeros términos
w=@(x,t) 0;
for n=1:10
w=@(x,t) w(x,t) + wi(x,t,n);
end
u=@(x,t) w(x,t)+v(x);
% Gráfica en 3D
close all
colormap jet
fsurf(u,[0 1 0 1], 'FaceColor','interp')
xlabel('x')
ylabel('t')
zlabel('u_1(x,t)')
Aunque a primera vista pueda parecer igual que la solución del problema original, si nos fijamos bien podemos observar que esta tarda más en alcanzar la solución estacionaria ya que esta curvada en un mayor intervalo, esto tiene sentido por tener un coeficiente de conductividad menor.
Ahora para apreciar mejor la diferencia vamos a representar ambas soluciones restándoles la estacionaria en x=1/2.
%Definimos las funciones
wi1=@(x,t,k) ((2*(-1)^(k))/(k*pi)).*sin(k.*pi.*x).*exp(-k^2.*pi^2.*t);
wi2=@(x,t,k) ((2*(-1)^(k))/(k*pi)).*sin(k.*pi.*x).*exp(-k^2.*pi^2.*t./2);
v=@(x) x;
%Sumamos los 10 primeros términos
w1=@(x,t) 0;
w2=@(x,t) 0;
for n=1:10
w1=@(x,t) w1(x,t) + wi1(x,t,n);
w2=@(x,t) w2(x,t) + wi2(x,t,n);
end
u1=@(x,t) w1(x,t)+v(x);
u2=@(x,t) w2(x,t)+v(x);
%Hacemos la diferencia con la solución estacionaria
dif1=@(t) u1(0.5,t)-v(0.5);
dif2=@(t) u2(0.5,t)-v(0.5);
% Gráfica
close all
t=linspace(0,1,100);
subplot(2,1,1)
plot(t,dif1(t),color='b',LineWidth=1.5)
xlabel('t')
ylabel('w(0.5,t)-v(0.5)')
title('k=1')
subplot(2,1,2)
plot(t,dif2(t),color='b',LineWidth=1.5)
xlabel('t')
ylabel('w(0.5,t)-v(0.5)')
title('k=1/2')
En estas gráficas podemos apreciar mejor que efectivamente esta solución tarda más en alcanzar la solución estacionaria.
4 Variación de las condiciones de frontera y las iniciales
En este apartado vamos a variar las condiciones de frontera y la inicial y a ver cómo varía la solución. El problema que vamos a estudiar ahora es el siguiente:
[math] \begin{cases} u_{2t}-u_{2xx}=0 & x\in(0,1),t\gt0\\ u_2(0,t)=0&t\gt0\\ u_2(1,t)=0&t\gt0\\ u_2(x,0)=max\{0,1-4|x-1/2|\}&x\in(0,1) \end{cases} [/math]
4.1 Resolución
Resolviendo este problema de la misma forma que los anteriores, obtenemos que la solución estacionaria [math] v_2(x)=0 [/math] y que la solución es: [math] u_2(x,t)=\sum_{k=1}^\infty c_k\sin(k\pi x)e^{-k^2\pi^2 t} [/math]. A diferencia de en los otros casos los [math]c_k[/math] los vamos a calcular aproximando las integrales por el método del trapecio
u0=@(x) max(0,1-4*abs(x-1/2));
uk=@(x,t,k) sin(k.*pi.*x).*exp(-k.^2*pi.^2.*t);
close all
a=0; b=1;
x=linspace(a,b,100);
t=linspace(a,b,100);
N=length(x)-1;
h=(b-a)/N;
w=ones(N+1,1);
w(1)=1/2; w(N+1)=1/2;
f_n=zeros(length(x),length(t));
for i=1:length(t)
for k=1:100
g=u0(x)*sin(k*pi*x)';
a_k=2*h*w'*g;
f_n(i,:)= f_n(i,:)+a_k.*uk(x,t(i),k);
end
end
colormap jet
surf(x,t,f_n,'FaceColor','interp')
xlabel('x')
ylabel('t')
zlabel('u_2(x,t)')
Al igual que en los casos anteriores, podemos ver que la solución alcanza la solución estacionaria.
4.2 Flujo en los extremos
Ahora, igual que para el problema original, vamos a ver cómo se comporta el flujo en este problema.
u0=@(x) max(0,1-4*abs(x-1/2));
uk=@(x,t,k) sin(k.*pi.*x).*exp(-k.^2*pi.^2.*t);
duk=@(x,t,k) -k*pi.*cos(k.*pi.*x).*exp(-k.^2*pi.^2.*t);
close all
a=0; b=1;
x=linspace(a,b,100);
t=linspace(a,b,100);
N=length(x)-1;
h=(b-a)/N;
w=ones(N+1,1);
w(1)=1/2; w(N+1)=1/2;
q0=@(t)0;
q1=@(t)0;
for k=1:100
g=u0(x)*sin(k*pi*x)';
a_k=2*h*w'*g;
q0=@(t) q0(t)+a_k.*duk(0,t,k);
q1=@(t) q1(t)+a_k.*duk(1,t,k);
end
subplot(2,1,1)
plot(t,q0(t),color='b',LineWidth=1.5)
xlabel('t')
ylabel('-u_{2x}(0,t)')
title('Flujo en el extremo izquierdo')
subplot(2,1,2)
plot(t,q1(t),color='b',LineWidth=1.5)
xlabel('t')
ylabel('-u_{2x}(1,t)')
title('Flujo en el extremo derecho')
4.3 Aislamiento de un extremo
[math] \begin{cases} u_{3t}-u_{3xx}=0 & x\in(0,1),t\gt0\\ u_3(0,t)=0&t\gt0\\ u_{3x}(1,t)=0&t\gt0\\ u_3(x,0)=max\{0,1-4|x-1/2|\}&x\in(0,1) \end{cases} [/math] [math] v_3(x,t)=0 [/math]
Solucion [math] u_3(x,t)=\sum_{k=1}^\infty c_k\sin((k+\frac{1}{2})\pi x)e^{-(k+\frac{1}{2})^2\pi^2 t} [/math]
5 Solución fundamental
La llamada solución fundamental de la ecuación del calor es la solución correspondiente a la condición inicial de un punto donde se concentra todo el calor.
- [math]\begin{cases} u_t - k u_{xx} = 0& (x, t) \in \mathbb{R} \times (0, \infty)\\ u(x,0)=\delta(x)& \end{cases}[/math]
donde [math]\delta[/math] es la función delta de Dirac, que vale 0 en todo punto menos en x, y su integral en toda la recta real es igual a 1.
La solución a este problema es la solución fundamental:
- [math]\Phi(x,t)=\frac{1}{\sqrt{4\pi kt}}\exp\left(-\frac{x^2}{4kt}\right).[/math]
%Parámetros
k = 1; % Coeficiente de conductividad
% Definición de los intervalos
x = linspace(-1, 1, 100);
t = linspace(10^-2, 1, 100);
% Crear la figura
figure;
% Inicialización de la matriz de soluciones
u = zeros(length(x), length(t));
% Cálculo de la solución fundamental
for i = 1:length(t)
u(:, i) = 1/sqrt(4 * pi * k * t(i)) * exp(-x.^2 / (4 * k * t(i)));
end
% Gráfica en 3D
colormap jet;
surf(x, t, u', 'FaceColor', 'interp');
title('Solución Fundamental de la Ecuación del Calor');
xlabel('x');
ylabel('t');
zlabel('Temperatura (u(x, t))');
ANIMACIÓN:
% Bucle para crear la animación
for i = 1:length(t)
u = 1/sqrt(4 * pi * k * t(i)) * exp(-x.^2 / (4 * k * t(i)));
% Actualizar la gráfica
plot(x, u, 'LineWidth', 2);
% Configurar el eje y la etiqueta
title(['Evolución de la Temperatura (t = ' num2str(t(i)) ')']);
xlabel('Posición (x)');
ylabel('Temperatura (u(x, t))');
% Establecer límites del eje y para mantener la escala
ylim([0, 3]);
% Pausa para controlar la velocidad de la animación
pause(0.08);
% Borrar la gráfica anterior para la siguiente iteración
if i < length(t)
clf;
end
end
De manera análoga, podemos extender la solución a dos dimensiones espaciales: [math] \Phi(x,y,t) = \frac{1}{4 \pi k t} \exp\left(-\frac{x^2 + y^2}{4 k t}\right).[/math]
% Parámetros
k = 1; % Difusividad térmica
% Definición de los intervalos en x y y
x = linspace(-1, 1, 100);
y = linspace(-1, 1, 100);
% Crear una malla 2D
[X, Y] = meshgrid(x, y);
% Tiempos para la animación
t_values = linspace(0.001, 0.1, 100);
% Crear la figura
figure;
% Bucle para crear la animación
for i = 1:length(t_values)
% Calcular la solución fundamental en 2D
u = 1/(4 * pi * k * t_values(i)) * exp(-(X.^2 + Y.^2) / (4 * k * t_values(i)));
% Actualizar la gráfica
surf(X, Y, u, 'EdgeColor', 'none');
% Configurar el eje y las etiquetas
title(['Solución Fundamental en 2D (t = ' num2str(t_values(i)) ')']);
xlabel('Posición x');
ylabel('Posición y');
zlabel('Temperatura u(x, y, t)');
% Establecer límites del eje y para mantener la escala
zlim([0, 20]);
% Pausa para controlar la velocidad de la animación
pause(0.2);
% Borrar la gráfica anterior para la siguiente iteración
if i < length(t_values)
clf;
end
end
5.1 Soluciones autosemejantes
Lo bueno que tiene la solución fundamental es que no hemos tenido que restringirla a una región acotada del espacio, sino que es solución para todo [math] x \in \mathbb{R} [/math]. No estaría de más preguntarse si existen otras soluciones con esta propiedad. La solución fundamental tiene la propiedad de ser autosemejante, es decir, su gráfica al variar el tiempo mantiene siempre la misma forma (una gaussiana).
Podemos buscar otras soluciones con esta propiedad considerando el cambio de variable [math]\xi=\frac{x}{\sqrt{t}}[/math]. Se comprueba que esta combinación de variables es invariante respecto de dilataciones parabólicas y además no tiene dimensión (PÁG 44 DEL LIBRO).
Imponemos que la función de una variable [math] U(\xi) [/math] cumpla la EDP, y llegamos a la EDO:
Cuya solución, deshaciendo el cambio de variable, es:
con [math] A,B [/math] constantes a determinar. Para ello, consideraremos las siguientes restricciones, que involucran una región espacial no acotada: [math] x\gt0 [/math]
[math] \left\{ \begin{aligned} &u(0, t) = 1 \\ &u(x, 0) = 0\\ &lim_{x \to \infty } u(x,t)=0 \end{aligned} \right. [/math]
Con las que obtenemos [math] A=1 [/math] y [math] B=0 [/math], y por tanto la solución a este problema:
% Definir el intervalo espacial
x = linspace(0, 5, 100);
% Definir el intervalo de tiempo
t_values = linspace(0.01, 1, 100);
% Crear la figura
figure;
% Bucle para crear la animación
for i = 1:length(t_values)
% Calcular la solución
u = exp(-x / (2 * sqrt(t_values(i))));
% Actualizar la gráfica
plot(x, u, 'LineWidth', 2);
% Configurar el eje y las etiquetas
title(['Solución de la Ecuación del Calor (t = ' num2str(t_values(i)) ')']);
xlabel('Posición (x)');
ylabel('Temperatura (u(x, t))');
% Establecer límites del eje y para mantener la escala
ylim([0, 1]);
% Pausa para controlar la velocidad de la animación
pause(0.1);
% Borrar la gráfica anterior para la siguiente iteración
if i < length(t_values)
clf;
end
end
5.2 La convolución
Otra forma de obtener soluciones de la ecuación del calor es usando la solución fundamental [math] \Phi(x,t) [/math] y aplicando una convolución.
- [math] u(x,t) = \int_{-\infty}^{\infty} \Phi(x-y,t) u_0(y) dy[/math]
Donde [math] u_0(x) [/math] es el dato inicial, es decir, [math] u(x,0) = u_0(x) [/math] para [math] x \in \mathbb{R} [/math]
Hemos visto antes que la solución fundamental representaba la respuesta del sistema a un impulso unitario (delta de Dirac) en un punto y en un instante de tiempo.
Intuitivamente, la convolución realiza una integración ponderada de las contribuciones locales del dato inicial en cada punto y momento, proporcionando así la evolución completa de la temperatura de acuerdo con la ecuación del calor.
Para ilustrarlo, vemos qué forma tiene la solución para el dato incial: [math] u_0(x) = 1_{[-1,1]} [/math]
(Falta cambiar el comando integral() por el de trapecio)
% Parámetros
k = 1; % Difusividad térmica
% Definición de los intervalos
x = linspace(-3, 3, 100);
t_values = linspace(0.0005, 0.3, 100);
% Función dato inicial
f = @(x) double(x >= -1 & x <= 1);
% Solución fundamental
G = @(x, t) 1/sqrt(4 * pi * k * t) * exp(-x.^2 / (4 * k * t));
% Crear la figura
figure;
% Bucle para crear la animación
for i = 1:length(t_values)
% Calcular la convolución
u = zeros(size(x));
for j = 1:length(x)
integrand = @(y) G(x(j) - y, t_values(i)) .* f(y);
u(j) = integral(integrand, -inf, inf);
end
% Actualizar la gráfica
plot(x, u, 'LineWidth', 2);
% Configurar el eje y la etiqueta
title(['Evolución de la Temperatura (t = ' num2str(t_values(i)) ')']);
xlabel('Posición (x)');
ylabel('Temperatura (u(x, t))');
% Establecer límites del eje y para mantener la escala
ylim([0, 1.1]);
% Pausa para controlar la velocidad de la animación
pause(0.01);
% Borrar la gráfica anterior para la siguiente iteración
if i < length(t_values)
clf;
end
end
