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. Para ello, vamos a ver cómo varía la solución en función de la variación de distintos parámetros como las condiciones iniciales o el coeficiente de conductividad. Además, veremos qué ocurre cuando suponemos que la longitud de la varilla no está acotada.
Contenido
1 Preliminares
Antes de empezar a desarrollar el problema necesitamos conocer algunos conceptos que nos van a permitir modelizarlo y sacar sus ecuaciones.
- Como el calor es una forma de energía, vamos a utilizar el principio de conservación de energía. Este establece que en un volumen de control [math]V[/math] arbitrario, la variación de energía a lo largo del tiempo se debe al balance neto del flujo de calor que atraviesa la frontera de [math]V[/math] ([math]\partial V[/math]) más una producción externa.
Donde [math] e(x,t) [/math] representa la densidad de energía calorífica y por tanto la primera integral se corresponde con la energía total que hay en [math]V[/math], [math] \overrightarrow q [/math] es el flujo de calor y como la segunda integral es una integral de superficie, representa el flujo que atraviesa la frontera de [math] V [/math], y [math] r [/math] es lo que aporta una fuente externa.
- También vamos a necesitar expresar el flujo y la energía en función de la temperatura. Para lo primero utilizamos la ley de Fourier, que 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, y u la temperatura. Y para lo segundo vamos a suponer que la energía es proporcional a la temperatura.
Siendo [math] c [/math] el calor específico.
- Denominamos coeficiente de difusión a [math] D=\frac{k}{c}. [/math]
- Principio del máximo: Sea [math] u\in C^{2,1}(Q_T)\cap C(\overline Q_T) [/math] tal que [math] u_t-\nabla u\leq 0 [/math] en [math] Q_T [/math]. Entonces [math] u [/math] alcanza su máximo en la frontera parabólica:
Análogamente, se cumple con el mínimo.
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 la temperatura se mantiene 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-Du_{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]
Siendo [math]u(x,t)[/math] función que representa la temperatura de la varilla a lo largo del tiempo. Vamos a suponer que [math]D=1[/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 sistema, 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 del problema original, 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 sistemas 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); %Solución del problema homogeneizado
v=@(x) x; %Solución estacionaria
%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); %Solución del problema
% Gráfica en 3D
close all
colormap jet
fsurf(u,[0 1 0 1], 'FaceColor','interp')
%Ejes
xlabel('x')
ylabel('t')
zlabel('u(x,t)')
En la gráfica podemos observar que la solución del problema alcanza en poco tiempo la solución estacionaria, lo cual tiene sentido ya que estamos utilizando un coeficiente de difusión relativamente alto
2.3 Flujo en los extremos
Ahora, para entender mejor cómo funciona la solución, vamos a ver cómo es su flujo. Para ello, vamos a graficar el flujo en los extremos teniendo en cuenta que, por la ley de Fourier, [math] q=-k u_x [/math].
%Definimos la función
dwi=@(x,t,k) (2*(-1)^(k+1)).*cos(k.*pi.*x).*exp(-k^2.*pi^2.*t); %Derivada de la solución del problema homogéneo
dv=@(x) 1; %Derivada de la solución estacionaria
%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); %Derivada de la solución en 0
du1=@(t) dw1(t)-dv(1); %Derivada de la solución en 1
%Definimos el vector de tiempo
t=linspace(0,1,100);
% Gráfica en x=0
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')
%Gráfica en x=1
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 extremo 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 la variación es más moderada, y por el contrario, en el extremo derecho es más drástica. 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]
Resolviéndolo de la misma forma que el anterior obtenemos que la estacionaria no varía y es [math] v_1(x)=x [/math], y que la solución es:
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); %Solución del problema homogéneo
v=@(x) x; %Solución estacionaria
%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); %Solución del problema
% Gráfica en 3D
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 está 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); %Solución del problema homogéneo con D=1
wi2=@(x,t,k) ((2*(-1)^(k))/(k*pi)).*sin(k.*pi.*x).*exp(-k^2.*pi^2.*t./2); %Solución del problema homogéneo con D=1/2
v=@(x) x; %Solución homogénea
%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); %Solución del problema con D=1
u2=@(x,t) w2(x,t)+v(x); %Solución del problema con D=1/2
%Hacemos la diferencia con la solución estacionaria en x=1/2
dif1=@(t) u1(0.5,t)-v(0.5);
dif2=@(t) u2(0.5,t)-v(0.5);
%Definimos el vector de tiempos
t=linspace(0,1,100);
% Gráfica de la solución del problema con D=1
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')
% Gráfica de la solución del problema con D=1/2
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. Suponemos que en el extremo derecho ahora la temperatura se mantiene a 0º, y que en el momento inicial la temperatura sigue la siguiente función: [math]max\{0,1-4|x-1/2|\}.[/math] Teniendo todo esto en cuenta, 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 es [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 que se utilizan para obtenerlos por el método del trapecio.
Haciendo todo esto, la gráfica de la solución queda de la siguiente forma:
%Definimos las funciones
u0=@(x) max(0,1-4*abs(x-1/2)); %Condición inicial
uk=@(x,t,k) sin(k.*pi.*x).*exp(-k.^2*pi.^2.*t); %Solución del problema
%Aplicamos la regla del trapecio
a=0; b=1; %Extremos del intervalo
x=linspace(a,b,100); %Vector de espacio
t=linspace(a,b,100); %Vector de tiempo
N=length(x)-1; %Número de puntos
h=(b-a)/N;
w=ones(N+1,1); %Vector de pesos
w(1)=1/2; w(N+1)=1/2;
f_n=zeros(length(x),length(t));
%Suma de los 10 primeros términos en todo el vector de tiempo
for i=1:length(t)
for k=1:10
g=u0(x)*sin(k*pi*x)';
a_k=2*h*w'*g; %Resultado de la regla del trapecio (c_k)
f_n(i,:)= f_n(i,:)+a_k.*uk(x,t(i),k);
end
end
%Gráfica en 3D
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 acaba alcanzando la solución estacionaria. Además, podemos observar que la solución siempre es positiva, lo que implica que la difusión transporta energía con velocidad infinita.
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. Para ello, vamos a representar el flujo en ambos extremos.
%Definimos las funciones
u0=@(x) max(0,1-4*abs(x-1/2)); %Condición inicial
uk=@(x,t,k) sin(k.*pi.*x).*exp(-k.^2*pi.^2.*t); %Solución del problema
duk=@(x,t,k) -k*pi.*cos(k.*pi.*x).*exp(-k.^2*pi.^2.*t); %Derivada de la solución
%Aplicamos el método del trapecio
a=0; b=1; %Extremos de los intervalos
x=linspace(a,b,100); %Vector de espacio
t=linspace(a,b,100); %Vector de tiempo
N=length(x)-1; %Número de puntos
h=(b-a)/N;
w=ones(N+1,1); %Vector de pesos
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; %Solución de la regla del trapecio (c_k)
q0=@(t) q0(t)+a_k.*duk(0,t,k); %Flujo en x=0
q1=@(t) q1(t)+a_k.*duk(1,t,k); %flujo en x=1
end
%Gráfica en x=0
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')
%Gráfica en x=1
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')
En ambos extremos el flujo empieza siendo cero, cosa que tiene sentido ya que tanto los extremos como los puntos de la barra cercanos a los extremos están a 0º. Posteriormente ambos flujos aumentan hacia fuera de la barra, esto es coherente ya que el calor en el centro de la barra en el momento inicial provoca que los laterales de la barra se vayan calentando y por tanto el flujo vaya hacia los extremos que se mantienen a 0º.
Una vez avanza el tiempo podemos observar que ambos flujo se estacionan en un valor, y esto coincide con el momento en el que la gráfica de la solución alcanza la solución estacionaria. Además, cabe resaltar que este valor es 0, lo cual tiene sentido ya que la solución estacionaria hace que toda la barra se quede a 0º.
4.3 Aislamiento de un extremo
En este apartado vamos a ver cómo se comporta la solución cuando aislamos térmicamente el extremo derecho, es decir, cuando el flujo en ese punto es cero. Teniendo esto en cuenta, el problema queda de la siguiente forma:
[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]
Resolviéndolo de la misma forma que en los casos anteriores obtenemos que la solución estacionaria es [math] v_3(x,t)=0 [/math] y que la solución del problema es [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]. Ahora, para calcular los coeficientes [math] c_k [/math], vamos a exigir que esta solución cumpla la condición inicial, es decir, [math] u_3(x,t)=\sum_{k=1}^\infty c_k\sin((k-\frac{1}{2})\pi x)= max\{0,1-4|x-1/2|\}[/math].
A diferencia de los otros casos, las autofunciones no las podemos expresar en función de la base trigonométrica, pero aun así forman un conjunto completo en [0,1] y vamos a comprobar que son ortogonales entre sí:
Como [math]\{ \sin((n-\frac{1}{2})\pi x) \}_{n\in\mathbb{N}}[/math] cumplen estas características, forman una base Hilbertiana y por tanto pueden aproximar cualquier función y los [math] c_k [/math] se pueden calcular igual que los coeficientes de Fourier pero utilizando esta nueva base. Para hacer esto vamos a utilizar la regla del trapecio.
Teniendo todo esto en cuenta, vamos a representar la función para entender mejor cómo funciona este sistema.
%Definimos las funciones
u0=@(x) max(0,1-4*abs(x-1/2)); %Condición inicial
uk=@(x,t,k) sin((k-0.5).*pi.*x).*exp(-(k-0.5).^2*pi.^2.*t); %Solución del problema
%Aplicamos la regla del trapecio
a=0; b=1; %Extremos de los intervalos
x=linspace(a,b,100); %Vector de espacio
t=linspace(a,b,100); %Vector de tiempo
N=length(x)-1; %Número de puntos
h=(b-a)/N;
w=ones(N+1,1); %Vector de nudos
w(1)=1/2; w(N+1)=1/2;
f_n=zeros(length(x),length(t));
%Sumamos los 100 primero términos en todo el vector t
for i=1:length(t)
for k=1:100
g=u0(x)*sin((k-0.5).*pi.*x)';
a_k=2*h*w'*g; %Solución de la regla del trapecio (c_k)
f_n(i,:)= f_n(i,:)+a_k.*uk(x,t(i),k);
end
end
%Gráfica en 3D
colormap jet
surf(x,t,f_n,'FaceColor','interp')
xlabel('x')
ylabel('t')
zlabel('u_3(x,t)')
Podemos observar que la gráfica es igual que la del problema anterior, excepto alrededor de x=1. Esto tiene sentido ya que lo único que hemos variado es la condición frontera en ese extremo. Además, podemos observar que la gráfica en esa zona es más recta, esto es debido a que la derivada es 0.
4.3.1 Flujo en los extremos
Ahora, igual que en todos los problemas anteriores, para entender mejor la solución vamos a analizar el flujo viendo cómo se comporta en los extremos mediante la siguiente gráfica.
%Definimos las funciones
u0=@(x) max(0,1-4*abs(x-1/2)); %Condición inicial
uk=@(x,t,k) sin((k-0.5).*pi.*x).*exp(-(k-0.5).^2*pi.^2.*t); %Solución del problema
duk=@(x,t,k) -(k-0.5)*pi.*cos((k-0.5).*pi.*x).*exp(-(k-0.5).^2*pi.^2.*t); %Derivada de la solución
%Aplicamos el método del trapecio
a=0; b=1; %Extremos de los intervalos
x=linspace(a,b,100); %Vector de espacio
t=linspace(a,b,100); %Vector de tiempo
N=length(x)-1; %Número de puntos
h=(b-a)/N;
w=ones(N+1,1); %Vector de pesos
w(1)=1/2; w(N+1)=1/2;
q0=@(t)0;
q1=@(t)0;
for k=1:100
g=u0(x)*sin((k-0.5)*pi*x)';
a_k=2*h*w'*g; %Solución de la regla del trapecio (c_k)
q0=@(t) q0(t)+a_k.*duk(0,t,k); %Flujo en x=0
q1=@(t) q1(t)+a_k.*duk(1,t,k); %flujo en x=1
end
%Gráfica en x=0
subplot(2,1,1)
plot(t,q0(t),color='b',LineWidth=1.5)
xlabel('t')
ylabel('-u_{3x}(0,t)')
title('Flujo en el extremo izquierdo')
%Gráfica en x=1
subplot(2,1,2)
plot(t,q1(t),color='b',LineWidth=1.5)
xlabel('t')
ylabel('-u_{3x}(1,t)')
title('Flujo en el extremo derecho')
En la gráfica del extremo izquierdo podemos observar que es igual que la anterior, al igual que pasa con la gráfica de la solución. Esto se debe a que la condición de contorno en ese extremo no la hemos variado.
Por el contrario, en este extremo la gráfica es distinta debido a que al aislar este extremo no va a haber flujo y por eso la gráfica es prácticamente 0 en todo momento.
5 Principio del máximo
Ahora, para entender mejor cómo funciona el principio del máximo, vamos a comprobar si este se cumple en las soluciones de los problemas.
En las dos primeras gráficas, podemos observar que el mínimo está en la frontera t=0 y el máximo en x=1 y en las otras dos podemos ver que el máximo está en t=0 y el mínimo en las fronteras de x. Por lo que efectivamente, en todas las soluciones que hemos calculado en este artículo se cumple el principio del máximo.
6 Solución fundamental
En este apartado vamos a estudiar la ecuación del calor cuando la variable espacial no está acotada. Por lo que el sistema queda de la siguiente forma:
- [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 sistema es la llamada solución fundamental de la ecuación del calor, correspondiente a la condición inicial de un punto donde se concentra todo el calor y su fórmula es:
- [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))');
De manera análoga, podemos extender esta 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; % Coeficiente de conductividad
% 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
colormap jet
surf(X, Y, u, 'EdgeColor', 'none');
% Configurar el eje y las etiquetas
title(['Solución Fundamental en 2D (t = ' num2str(t_values(i)) ')']);
xlabel('x');
ylabel('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
6.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ágina 44 Partial Differential Equations in Action. S. Salsa, G. Verzini).
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
6.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]
% Definición de los intervalos
x = linspace(-5, 5, 1000);
t_values = linspace(0.01, 1, 100);
% Función dato inicial
f = @(x) double(x >= -1 & x <= 1);
% Solución fundamental
G = @(x, t) 1/sqrt(4 * pi * t) * exp(-x.^2 / (4 * t));
% Crear la figura
figure;
% Bucle para crear la animación
for i = 1:length(t_values)
% Calcular la convolución utilizando el método del trapecio
u = zeros(size(x));
for j = 1:length(x)
integrand = @(y) G(x(j) - y, t_values(i)) .* f(y);
u(j) = trapz(x, integrand(x));
end
% Actualizar la gráfica
plot(x, u, 'LineWidth', 2);
% Configurar el eje y las etiquetas
title(['Convolución con Trapecio (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.1);
% Borrar la gráfica anterior para la siguiente iteración
if i < length(t_values)
clf;
end
end
