Ecuación del calor ( ALA )

Trabajo realizado por estudiantes
Título	Ecuación del Calor.
Asignatura	EDP
Curso	2023-24
Autores	Andrea Navarro, Lucía Amores, Aitana Guill.
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

En este documento se pretende mostrar al lector como la ecuación del calor en una dimensión describe el fujo de calor [math] u(x,t) [/math] Para ello estudiaremos distintas condiciones frontera e iniciales en una barra metálica que ocupa un intervalo [0,1].

2 Sistema no homogéneo

En este primer caso nos centraremos en una barra metálica que comienza estando a 0 °C y cuyas temperaturas al principio y al final de son dos constantes distintas. En concreto, consideraremos que la temperatura en la posición x = 0 es nula, y sin embargo, en x = 1 la sube un grado. Asimismo, estudiaremos la ecuación del calor cuya conductividad térmica, k, y calor específico consideraremos 1. Todo esto se traduce en el sistema no homogéneo,

[math] \left\{ \begin{array}{ll} u_{t}(x,t) -u_{xx}(x,t)=0, \hspace{5mm} x\in[0,1] \hspace{3mm} t\gt0 \\ u(x,0)=0, \hspace{5mm} x\in[0,1] \\ u(0,t)=0,\hspace{3mm} t\gt0\\ u(1,t)=1,\hspace{3mm} t\gt0 \end{array} \right. [/math]

2.1 Solución Estacionaria

La resolución del sistema anterior se basa en el método de separación de variables perteneciente a la teoría de resolución ecuaciones diferenciales, lo cual carece de interés en este documento. Por limpieza en la lectura, este procedimiento no se incluirá, si embargo, hay un paso previo al método que cabe incluir. La resolución por separación de variables requiere que el sistema sea homogéneo. Esta modificación en el sistema original la conseguimos haciendo uso de la que se conoce como solución estacionaria. Esta se alcanza cuando ha pasado un tiempo infinito ([math] t \rightarrow \infty [/math]) y considerando por tanto, que el flujo del calor ha dejado de depender del tiempo, [math] u(t,x) \sim v(x) [/math]. Haciendo los respectivos cálculos es fácil llegar a que la solución estacionaria es [math] v(x)=x [/math] con [math] x\in[0,1] [/math]. La cual gráficamente muestra una bajada de temperatura en el espacio y tiempo.

Las gráficas de la derecha han sido dibujadas en Matlab implementando el siguiente código.

Solución estacionaria

clear
close all
clc

% La EDO de la EDP estacionaria viene definida por:
% -v_xx=0 con condiciones frontera: v(0)=0 y v(1)=1
% La solución viene dada por:
v=@(x) x;

x=0:0.001:1; %longitud

figure
subplot(2,1,1)
plot(x,v(x),'b',LineWidth=1.5)
title('Solución estacionaria del sistema de EDP en ℝ^2')
legend('v(x) = x', location='northwest')
xlabel('x')
subplot(2,1,2)
[X,T]=meshgrid(x,0:0.001:1);
V=v(X);
mesh(X,T,V)
xlabel('x')
ylabel('t')
title('Solución estacionaria del sistema de EDP en ℝ^3')

2.2 Solución del sistema no homogéneo

Una vez hallada la solución estacionaria, el método de separación de variables nos devuelve el candidato a solución, [math] u_k(x,t)=x-c_k\cdot e^{-k^2\pi^2 t}\cdot \sin{(k\pi x)}[/math] Donde [math] c_{k} = \frac{2(-1)^k}{k\pi} [/math] son los coeficientes de Fourier hallados mediante aproximación impar. Al igual que con la solución estacionaria, podemos ver que la temperatura va decreciendo en el espacio y tiempo.

Solución del sistema no homogéneo

Más aún, es fácil observar como la tras el paso del tiempo, esta solución se va aproximando a la estacionaria.

La gráfica ha sido dibujada con el siguiente código de Matlab.

clear
close all
clc

% Creamos un vector con los intervalos a tratar:
x=0:0.001:1; %longitud de la varilla
t=0:0.001:1; %tiempo

% Los coeficientes de Fourier vienen dados por:
K=10;
c=zeros(K,1);
for k=1:K
    c(k)=2*(-1)^(k+1)/(k*pi);
end

% La solución del sistema de EDP viene dada por:
w=@(x,t)0;
for k=1:K
    w=@(x,t) w(x,t) + c(k).*exp(-k^2*pi^2.*t).*sin(k.*pi.*x);
end

% La solución estacionaria:
v=@(x)x;

% Definimos una malla de puntos:
[X,T]=meshgrid(x,t);

% Valores de w en forma de malla de puntos:
W=w(X,T);


% Valores de v en la malla de puntos:
V=v(X);

% Valores de U en la malla de puntos:
U=V-W;

% Representamos gráficamente la solución de u(x,t):
figure
mesh(X,T,U)
xlabel('x')
ylabel('t')
zlabel('temperatura')
title('Solución de u(x,t)')

2.3 Flujo del calor entrante y saliente

Una vez vista tanto la solución del sistema como el estado estacionario de la barra, en este apartado pretendemos entender de formas más intuitiva como fluye el flujo de calor por la barra metálica. Para ello es especialmente importante tener en cuenta los lados izquierdo y derecho de la barra, ya que es a través de ellos por donde sale y entra el flujo respectivamente. Introduciendo el siguiente código en Matlab conseguimos una gráfica que nos ayudará a estudiar estos movimientos.

clear
close all
clc

% Intervalos a tratar:
t=0:0.001:1; %tiempo

% Coeficientes de Fourier:
K=10;
c=zeros(K,1);
for k=1:K
    c(k)=2*(-1)^(k+1)/(k*pi);
end

% Derivada de w(x,t) respecto de x:
dw_x=@(x,t)0;
for k=1:K
    dw_x=@(x,t)dw_x(x,t)+c(k)*exp(-k^2*pi^2*t).*(k*pi).*cos(k*pi.*x);
end

% Derivada de v(x,t) respesto de x:
Dv=@(x)1;

% El flujo viene dado por -ku_x, siendo k el coeficiente de conductividad
% del calor y u_x la derivada de la solución u(x,t) respecto de x:

kk=1; %coeficiente de conductividad del calor

% El flujo entrante x=0:
DU0=Dv(0)-dw_x(0,t);   %derivada de u en x=0
flujo_x_0 = -kk.*DU0; %flujo entrante

% El flujo saliente x=1:
DU1=Dv(1)-dw_x(1,t); % derivada de u en x=1
flujo_x_1=-kk.*DU1; %flujo saliente

% Representación gráfica de los flujos:
figure
subplot(1,3,1)
plot(t,flujo_x_0,'b','LineWidth',1.5)
title('Flujo Entrante')
grid on
grid minor
axis square
subplot(1,3,2)
plot(t,flujo_x_1,'r','LineWidth',1.5)
title('Flujo Saliente')
grid on
grid minor
axis square
subplot(1,3,3)
hold on
plot(t,flujo_x_0,'b--','LineWidth',1.5)
plot(t,flujo_x_1,'r--','LineWidth',1.5)
hold off
grid on
grid minor
axis square
title('Flujos')
legend('Flujo entrante','Flujo saliente',Location='southeast')

Flujos entrante y saliente

A simple vista podemos ver que al principio el flujo entrante tiene una variación mucho menos marcada que la del flujo saliente. Esto simplemente se debe a que en x=0 la temperatura ambiente es la misma que la de la barra metálica, mientras que en x=1 la barra está a un grado más. Además, cabe mencionar que ambos flujos se estabilizan cuando el flujo llega a -1 porque es aquí cuando la solución del sistema alcanza el estado estacionario.

Si nos centramos en los valores que toma el flujo saliente, podemos apreciar que este tiende a hacerse más negativo. Esto físicamente se debe a que el flujo tiende a ir hacia donde hay menos temperatura, y el extremo izquierdo se mantiene a cero grados mientras que el resto de la barra se calienta. Por el otro lado de la barra tenemos que en los primeros instantes el flujo cada vez toma valores más grandes, entrando así menos flujo por x=1. Sin embargo, una vez el lado izquierdo se va calentando, comienza a entrar más flujo hasta estabilizarse por completo.

2.4 ¿Qué pasa si la tomamos [math] k = \frac{1}{2} [/math]?

De momento en la ecuación del calor hemos considerado todas las constantes como 1 pero, ¿qué pasaría si la conductividad térmica disminuye?, ¿podremos apreciar algún cambio grande en la solución final?. Para ello tomamos [math] k = \frac{1}{2} [/math] y resolvemos el mismo sistema pero esta vez con ecuación,

[math] u_{t}(x,t) - \frac{1}{2}u_{xx}(x,t)=0, \hspace{5mm} x\in[0,1] \hspace{3mm} t\gt0 . [/math]

Homogeneizar dicho sistema nos devuelve la solución estacionaria [math] v(x)_{\frac{1}{2}} = x[/math], que gráficamente tiene la forma,

Solución estacionaria

De la misma forma que en el sistema con [math] k = 1 [/math], obtenemos la solución final,

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

siendo [math]c_k=\cfrac{2\cdot(-1)^{k+1}}{k\pi}[/math].

La cual gráficamente se comporta de la siguiente manera

Solución de u(x,t) con conductividad 1/2

Comparando esta gráfica con la correspondiente de [math] u(x,t) [/math] no somos capaces a simple vista de encontrar ninguna diferencia.

Una forma más visual de apreciar el efecto de la disminución de [math]k[/math] es comparando las soluciones de los sistemas con sus estacionarias correspondientes. Es decir, vamos a ver con cuanta "rapidez" cada sistema alcanzará su correspondiente solución estacionaria. Esto los haremos calculando las diferencias entre la solución y el estado estacionario para [math] x =\frac{1}{2} [/math] en ambos casos; [math] u(\frac{1}{2},t) = v(\frac{1}{2}) [/math] y [math] u_{\frac{1}{2}}(\frac{1}{2},t) = v_{\frac{1}{2}}(\frac{1}{2}) [/math].

Todas las gráficas se han representado con el siguiente código de Matlab:

derecha

clear
close all
clc

% Intervalos a tratar:
x=0:0.001:1; %longitud
t=0:0.001:1; %tiempo

% Número que nos indica los K primeros términos de la serie:
K=10;

% Solución estacionaria:
v=@(x)x;

% Calculamos los coeficientes de Fourier del sistema de EDP (homogeneo):
c=zeros(K,1);
for k=1:K
    c(k)=2*(-1)^(k+1)/(k*pi);
end

% Solución calculada del sistema de EDP (homogeneizado):
w=@(x,t)0;
for k=1:K
    w=@(x,t) w(x,t)+ c(k).*exp(-k.^2.*pi^2.*t*1/2).*sin(k*pi.*x);
end

u=@(x,t) v(x)-w(x,t); %Solición con conductividad 1/2

% Representamos una malla de puntos:
[X,T]=meshgrid(x,t);    
W=w(X,T);
U=u(X,T);

% Representación gráfica de la solución estacionaria:
figure
subplot(2,1,1)
plot(x,v(x),'b',LineWidth=1.5)
xlabel('x');
title('Solución estacionaria del sistema de EDP en ℝ^2')
legend('v(x) = x', location='northwest')

subplot(2,1,2)
V=v(X);
mesh(X,T,V);
title('Solución estacionaria del sistema de EDP en ℝ^3')
xlabel('x')
ylabel('t')

% Representamos gráficamente:
figure
mesh(X,T,U)
xlabel('x')
ylabel('t')
zlabel('u(x,t)')
title('Solución de u_1_/_2 (x,t)')

% Diferencia entre la solución y el estado estacionario:

% Primero volvemos a calcular la solución para conductividad 1:
% Los coeficientes de Fourier vienen dados por:
cc=zeros(K,1);
for k=1:K
    cc(k)=2*(-1)^(k+1)/(k*pi);
end

% La solución del sistema de EDP viene dada por:
ww=@(x,t)0;
for k=1:K
    ww=@(x,t) ww(x,t) + cc(k).*exp(-k^2*pi^2.*t).*sin(k.*pi.*x);
end

uu=@(x,t) v(x)-ww(x,t);

figure
xx=1/2;
hold on
plot(t,v(1/2)-u(xx,t),'b--')
plot(t,v(1/2)-uu(xx,t),'r--')
grid on
legend('u_1_/_2(1/2,t) - u_e_s_t(1/2)','u_1(1/2,t) - u_e_s_t(1/2)')
xlabel('t')
ylabel('u_k(1/2,t) - u_e_s_t(1/2)')
title('Diferencias entre la solución t el estado estacionario en x=1/2')

Habiendo entendido el papel que juega la solución estacionaria resulta bastante intuitivo que esta se alcance en más tiempo en aquella barra con menor conductividad térmica. Esto es debido a que el calor pasará con mayor dificultad por la barra, ralentizando así la llegada a la solución estacionaria.

3 Sistema con cambios en la condición inicial

En este apartado vamos a estudiar el papel que juega la condición inicial imponiendo una condición que varíe en el espacio. En concreto vamos a considerar la condición inicial [math] u(x,0) = \max\{0, 1-4·|x-1/2|\} - x [/math]. Para facilitar el estudio vamos a considerar directamente un sistema homogéneo, es decir, [math] u(0,t) = u(1,t) = 0 [/math]. El sistema completo sería,

[math] \left\{ \begin{array}{ll} u_{t}(x,t) - u_{xx}(x,t)=0, \hspace{5mm} x\in[0,1] \hspace{3mm} t\gt0 \\ u(x,0)= \max\{0, 1-4·|x-1/2|\} - x, \hspace{5mm} x\in[0,1] \\ u(0,t)=0,\hspace{3mm} t\gt0\\ u_t(1,t)=0,\hspace{3mm} t\gt0 \end{array} \right. [/math]

Debido a que el sistema ya es homogéneo de partida, la solución estacionaria pasaría a ser [math] v(x) = 0 [/math] perdiendo así interés.

3.1 Solución del sistema

Tras aplicar el correspondiente método de separación de variables obtenemos la solución del sistema,

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

Donde los coeficientes de Fourier [math]c_k = 2\int_{0}^{1}\sin(k \pi x)u(x,0)[/math], se han calculado en Matlab aproximando las integrales por el método del trapecio,

Solución estacionaria

Solución de u(x,t)

clear
close all
format long
clc

% Intervalos con los que trabajaremos:
x=0:0.001:1; %longitud
t=0:0.001:1; %tiempo

% Cálculo de los coeficientes de Fourier:
K=10; % K primeros términos de la serie
f=@(x) max(0, 1-4*abs(x-1/2)); %condicion inicial
% Fórmula  del trapecio para aproximación:
c= zeros(K,1); %matriz para almacenar los coeficientes de Fourier
for k=1:K
    c(k)=2*trapz(x,f(x).*sin(k*pi*x));
end

% La solución del sistema de EDP viene dada por:
u=@(x,t)0;
for k=1:K
    u=@(x,t) u(x,t)+c(k).*exp(-k^2*pi^2*t).*sin(k*pi*x);
end

% Creamos una malla de puntos:
[X,T]=meshgrid(x,t);
U=u(X,T);

% Representamos gráficamente:
figure
mesh(X,T,U)
xlabel('x')
ylabel('t')
zlabel('u(x,t)')
title('Solución de u(x,t)')

Veamos a medida que pasa el tiempo nuestra solución se aproxima al estado estacionario [math]v(x) = 0[/math]. Otra cosa interesante a mencionar de esta gráfica es que la solución es estrictamente positiva si [math]t \gt 0[/math]. Esto físicamente se puede interpretar pensando que la difusión lleva una velocidad infinita de propagación.

3.2 Flujo del calor entrante y saliente

Al igual que en el sistema anterior, veamos como fluye el flujo de calor a lo largo del tiempo.

Flujos entrante y saliente

INTRODUCIR CODIGO En el extremo x=0 tenemos un gráfica muy similar a la descrita en el sistema anterior ya que esta condición inicial no ha cambiado. Por otro lado, la gráfica del flujo entrante se vuelve cero rápidamente debido a que el extremo x=1 está aislado por las condiciones impuestas.

4 Sistema con cambios en las condiciones frontera

En esta sección se llevará a cabo un estudio del comportamiento de la solución tras aislar térmicamente el extremo derecho de la varilla metálica. Esto implica un flujo nulo en ese punto, obteniendo el siguiente sistema:

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

4.1 Solución del sistema.

Teniendo que las condiciones frontera son nulas, se tiene que la solución estacionaria [math]v(x)=0 [/math]. Por lo que la solución del problema es:

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

siendo [math]c_k[/math] los coeficientes de Fourier calculados por aproximación a través del método del trapecio con Matlab. Finalmente, teniendo en cuenta todo se representa gráficamente la solución a continuación para estudiar su comportamiento:

Solución para los 10 primeros términos de la serie

Solución para los 100 primeros términos de la serie

clear
close all
format long
clc

% Intervalos con los que trabajamos:
x=0:0.001:1; %longitud
t=0:0.001:1; %tiempo

K=10; % K primeros términos de la serie

f=@(x)max(0,1-4.*abs(x-1/2)); %condición inicial

% Calculamos los coeficientes de Fourier mediante la fórmula del trapecio:
c=zeros(K,1);

for k=1:K
    c(k)=2*trapz(x,f(x).*sin((((2.*k-1).*pi)/2).*x));
end

% Solución de nuestro sistema de EDP:
u=@(x,t)0;
for k=1:K
    u=@(x,t) u(x,t) + c(k)*exp(-(2.*k-1)^2.*pi^2*t/4).*sin((2.*k-1).*pi/2.*x);
end

% Creamos una malla de puntos:
[X,T]=meshgrid(x,t);
U=u(X,T);

% Representamos gráficamente la solución de la EDP:
figure
mesh(X,T,U)
xlabel('x')
ylabel('t')
zlabel('u(x,t)')
title('Solución de u(x,t)')

El comportamiento de la solución es muy similar al del apartado anterior, en el que se ven modificadas las condiciones frontera e inicial; a diferencia del entorno de [math]x=1[/math], ya que es el punto en el que se ha visto modificada la varilla.

4.2 Flujos entrante y saliente

Se analizarán los flujos en los extremos de la varilla de forma análoga a las secciones anteriores.

Flujos.

El flujo entrante, es decir, el asociado al extremo izquierdo de la varilla es idéntico al estudiado en la sección anterior, donde se cambiaron las condiciones frontera e inicial. Mientras que el asociado al extremo derecho es prácticamente constante y muy cercano a 0, esto se debe a que la varilla se encuentra aislada térmicamente en el extremo derecho, por lo que no habrá flujo.

Las gráficas y los flujos se han calculado mediante el siguiente código de Matlab:

clear
close all
format long
clc

% Intervalos con los que trabajamos:
x=0:0.001:1; %longitud
t=0:0.001:1; %tiempo

K=10; % K primeros términos de la serie

f=@(x)max(0,1-4.*abs(x-1/2)); %condición inicial del sistema

% Calculamos los coeficientes de Fourier:
c=zeros(K,1);

for k=1:K
    c(k)=2*trapz(x,f(x).*sin((((2.*k-1).*pi)/2).*x));
end

% Creamos una malla de puntos:
[X,T]=meshgrid(x,t);
U=u(X,T);

% El flujo de u(x,t) es -k*u_x(x,t) siendo k el coeficiente de
% conductividad del calor, y u_x la derivada de u respecto de x:

% Calculamos la derivada de u:
du=@(x,t)0;
for k=1:K
    du=@(x,t) du(x,t) + c(k)*(2.*k-1).*pi/2*exp(-(2.*k-1)^2.*pi^2*t/4).*cos((2.*k-1).*pi/2.*x);
end

kk=1; %coeficiente de conductividad del calor

% El flujo entrante x=0:
DU0=du(0,t);   %derivada de u en x=0
flujo_x_0 = -kk.*DU0; %flujo entrante

% El flujo saliente x=1:
DU1=du(1,t); % derivada de u en x=1
flujo_x_1=-kk.*DU1; %flujo saliente

% Representación gráfica de los flujos:
figure
subplot(1,3,1)
plot(t,flujo_x_0,'b','LineWidth',1.5)
title('Flujo Entrante')
grid on
grid minor
axis square
subplot(1,3,2)
plot(t,flujo_x_1,'r','LineWidth',1.5)
title('Flujo Saliente')
grid on
grid minor
axis square
subplot(1,3,3)
hold on
plot(t,flujo_x_0,'b--','LineWidth',1.5)
plot(t,flujo_x_1,'r--','LineWidth',1.5)
hold off
grid on
grid minor
axis square
title('Flujos')
legend('Flujo entrante','Flujo saliente',Location='southeast')

5 Principio del máximo

A continuación se muestran las soluciones de todos los sistemas planteados en las secciones anteriores, con la finalidad de ver si se cumple el principio del máximo:

Soluciones de los sistemas anteriores.

Como podemos apreciar en las gráficas, todas las soluciones de los sistemas planteados cumplen el principio del máximo, alcanzando tanto su máximo como su mínimo en la frontera.

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

Anteriormente se han representado diferentes soluciones de la ecuación del calor en una dimensión. En este apartado, se profundizará en la existencia de una función matemática que describe como se propaga y disipa el calor en un medio a lo largo del tiempo, denominada solución fundamental de la ecuación del calor. En este primer caso, se realizará para dimensión 1, la cual toma la forma de una función Gaussiana. Sin embargo, posteriormente se ampliará esta a una segunda dimensión.

La expresión de esta es la siguiente:

[math] u(x,t) = \frac{1}{ \sqrt{4\pi t}}e^{-\frac{1}{4t}x^2} [/math]

Para su correcta comprensión, se va a realizar a continuación una representación de la solución fundamental de la ecuación del calor en dimensión 1, particularizando el dominio espacial en el intervalo [-1,1] y el dominio temporal en el intervalo [0.01,1]. Para ello, se ha creado el siguiente código, explicado en las anotaciones de él.

clc
close all
format long
f = @(x,t) exp(-(x.^2)./(4*t))./(sqrt(4*pi*t)); %Se define la función de la solución fundamental de la ecuación del calor.
%Se grafica la representación de esta solución
[X,T]=meshgrid(-1:0.001:1, 10^-2:0.001:1);
mesh(X,T,f(X,T));
xlabel('x');
ylabel('t');
zlabel('u(x,t)');
title('Solución fundamental de la ecuación del calor en dimensión 1');

La representación obtenida tras ejecutar el código anterior es la siguiente:

centro

Como se observa en la imagen, la gráfica alcanza un máximo en el intervalo espacial y temporal [-0.5,0.5]x[0.01,0.25], el cual representa una propagación superior del calor cuando las variables alcanzan dichos valores. Sin embargo, se observa que en el resto de nuestro dominio, la propagación del calor sufre cambios inferiores, generando gráficamente una superficie ligeramente plana.

7 Solución de la ecuación del calor mediante la convolución

A continuación se presenta una nueva manera de encontrar la solución de la ecuación del calor u_t-u_xx=0 en valores de x reales. Para ello, se explicará mediante el siguiente ejemplo.

Sea la condición inicial u₀(x)=1_[-1,1], la solución a la ecuación anterior viene dada por la convolución siguiente, donde K(x,t) es la solución fundamental explicada en apartados anteriores:

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

Este enunciado se asemeja al siguiente teorema propuesto:

7.1 Teorema

Sea f una función continua y acotada en todos los reales, se cumple que [math] u(x,t) = \int_{-\infty}^{\infty} K(x-y,t)·f(y) dy [/math] es una solución de la ecuación del calor u_t-u_xx=0 con condición incial u(x,0)=f(x), siendo K la solución fundamental del calor vista anteriormente.

Se puede observar que nuestro enunciado y el teorema presentan diferencias. En el primero no se cumple la condición de que nuestra función de la condición inicial sea continua en todos los reales. Aún así, se ha representado la solución obtenida, llegando a la conclusión de que en este caso se puede obtener la solución mediante la convolución.

Para ello, se ha realizado el siguiente código en el cual se evalúa y representa gráficamente la convolución para valores del tiempo iguales a t=0.001, t=0.01 y t=0.1. A pesar de que el dominio de la x sean todos los reales, en la representación se ha decidido realizar la gráfica en el intervalo [-3,3], pues con intervalos superiores se dificulta su observación:

clear
close all
format long
% Se define la solución fundamental de la ecuación del calor en dimensión 1
K=@(x,t) 1./sqrt(4*pi*t)*exp(-x.^2./(4*t));

% Se define la condición inicial en el intervalo dado
u0=@(x)(x>=-1 & x<=1);


% Se define el vector que almacena los instantes de tiempo en lo que queremos evaluar la solución.
t_valores = [0.001, 0.01, 0.1];

% Se define el vector que alamcena los 1000 puntos equidistantes en el dominio de la variable espacial. Para que la representación se observe de manera óptima se toma el siguiente dominio. Aún así se puede representa en toda la recta real
x = linspace(-3, 3, 1000); 
dx = x(2) - x(1);

% Cálculo de la solución en los diferentes valores del vector de tiempo
figure;
for i = 1:length(t_valores)
    t = t_valores(i);
    
    % Aproximación de la integral utilizando la fórmula del trapecio
    integral_approx = zeros(size(x));
    for j = 1:length(x)
        y = x(j);
        integral_approx(j) = sum(K(x - y, t) .* u0(x)) * dx; % Aplicación de la fórmula del trapecio
    end
    
    % Representación gráfica
    subplot(1, length(t_valores), i);
    plot(x, integral_approx);
    xlabel('x');
    ylabel('u(x, t)');
    title(['t = ', num2str(t)]);
end

Una vez ejecutado el código anterior, se obtiene una representación gráfica de la solución en cada uno de los instantes. A continuación se presenta y posteriormente se analiza.

centro

Analizando la representación gráfica anterior se observa que conforme el valor del tiempo decrece aproximándose a cero, la solución se parece cada vez más a la función de la condición inicial, la cual se alcanza cuando el valor del tiempo es cero.

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

En los apartado anteriores, se ha trabajado con la solución fundamental de la ecuación del calor de dimensión uno. Sin embargo, también se puede trabajar con esta solución en dimensiones superiores. En el caso de la dimensión 2, la solución fundamental de la ecuación del calor es la siguiente:

[math] u(x,t) = \frac{1}{ 4\pi t}e^{-\frac{x^2+y^2}{4t}} [/math]

Para una óptima comprensión de como se comporta se va a representar gráficamente para diversos instantes de tiempo en el intervalo espacial (x1, x2) = [−1, 1]². Para ello se ha implementado el siguiente código y se ha obtenido la gráfica que aparece a continuación.

%Se define el vector de 1000 puntos de la variable espacial
x = linspace(-1, 1, 1000);
[X1, X2] = meshgrid(x, x);

% Se crea el vector de los tiempos que se van a evaluar
t_valores = [0.001, 0.01, 0.1];

% Graficar la solución en diferentes instantes de tiempo
figure;
for i = 1:length(t_valores)
    t = t_valores(i);
    
    % Calcular la solución fundamental en dimensión 2
    u = 1./(4*pi*t) * exp(-(X1.^2 + X2.^2)./(4*t));
    
    % Graficar la solución en el instante de tiempo actual
    subplot(1, length(t_valores), i);
    
    mesh(X1, X2, u);
    xlabel('x1');
    ylabel('x2');
    axis square
   
    title(['t = ', num2str(t)]);
    
    colormap(jet); % Colormap para visualizar mejor
end

centro

Como se observa en la imagen, la solución fundamental de la ecuación del calor presenta un máximo, lo cual permite ver que el calor en un instante espacial se propaga de manera superior al resto. A su vez, se observa que cuanto más próximo sea el valor de t a cero, la solución es más singular.

9 Soluciones autosemejantes

A continuación, se realizará un breve análisis de otro tipo de soluciones denominadas autosemejantes. Estas son invariantes bajo ciertos cambios de escala en espacio y tiempo. Para comprenderlo de manerá óptima se va a generar un ejemplo. Una de estas soluciones es de la forma [math]u(x,t)=U(\frac{x}{\sqrt{t}})[/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]

Si se añaden las anteriores condiciones iniciales, se llega a la conclusión de que mediante una serie de cambios de variables la solución es la siguiente:

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

Se ha realizado el siguiente código para obtener la representación de esta solución:

u = @(x, t) -1 * erf(x ./ (2 .* sqrt(t))) + 1;


[x,t] = meshgrid(0:.01:5);
mesh(x,t,u(x,t))
xlabel('x');
ylabel('t');

centro

Tal y como se observa en la imagen, la temperatura decrece conforme la variable espacial aumenta y la temporal decrece. Se analiza como la temperatura tiende a cero mientras el tiempo se aproxima a este valor.