Ecuación del calor. Yan, Otelo, Mika
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Ecuación del calor. Grupo 6-A |
| Asignatura | EDP |
| Curso | 2023-24 |
| Autores | Miguel Cazorla Pedraza
Otelo Gallego Ayala Yan Wang |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
Contenido
1 Introducción
En este articulo primero estudiaremos la ecuación del calor en dominio acotado, y más tarde trataremos su solución fundamental en todo [math] \mathbb{R} [/math]. Además, graficaremos los resultados para entender mejor la teoría de estas ecuaciones y compararemos diversos casos.
2 Ecuación del calor I
Vamos a dibujar diferentes soluciones de la ecuación del calor en una dimensión. Para ello, partimos del siguiente enunciado:
[math] \textit{Se considera una varilla metálica que ocupa el intervalo [0,1] y que se encuentra aislada por su superficie lateral, de manera que la conducción de calor sólo se produce}[/math] [math] \textit{en la dirección longitudinal. La temperatura inicial de la varilla es } [/math] [math] 0^oC [/math]. [math] \textit{ En el extremo izquierdo se consigue mantener la temperatura a } 0^oC \textit{mientras que en el} [/math] [math] \textit{ derecho la temperatura es de } 1^oC [/math].
Suponiendo que tanto la conductividad térmica como el calor específico tienen un valor constante [math] 1 [/math], y llamando [math] u(x,t) [/math] a la función temperatura en función del espacio y el tiempo, tenemos que el anterior enunciado se ve modelizado por el siguiente sistema EDP:
Para dar solución a este sistema, resolveremos primero el problema resultante de tomar [math] t \rightarrow \infty [/math], que implica [math] \frac{du}{dt}(x,t) = 0 [/math] por lo que se denomina caso estacionario.
La solución estacionaria es [math] V(x) = x [/math]. Veamos cómo es su representación gráfica en 3D:
%Resumen: creamos una matriz que representa una malla de puntos (tiempo,
%espacio) y en cada una de las columnas introducimos el valor de x para
%todos los tiempos
%matriz
evaluaciones = zeros(100,100);
%límite temporal y vectores
T = 1;
t = linspace(0,T,100);
x = linspace(0,1,100);
%rellenamos la matriz
for j=1:100
evaluaciones(:,j) = x(j).*ones(100,1);
end
%representación gráfica
surf(t,x,evaluaciones')
title('Solución estacionaria')
xlabel('Tiempo')
ylabel('Espacio')
zlabel('Temperatura (°C)')
La solución del sistema original por separación de las variables que obtenemos es:
Y su representación gráfica:
%Resumen: la idea del programa es crear una matriz de pares (tiempo,
%espacio) y en cada uno de esos puntos calcular la aproximación de la
%solución mediante un bucle for que realiza la suma de términos de la serie
%creamos la malla en la que introduciremos las evaluaciones
%la condición frontera del extremo izquierdo queda impuesta por ser una
%matriz de ceros
evaluaciones = zeros(100,100);
%imponemos la condición frontera del extremo derecho
evaluaciones(:,100) = ones(100,1);
%limite temporal y vectores de tiempos y espacios
T = 1;
t = linspace(0,T,100);
x = linspace(0,1,100);
%coeficiente de conductividad térmica
cct = 1;
%en cada punto de la malla sumamos los diez primeros términos de la serie
%obtenida como solución. Después, sumamos también la solución estacionaria,
%que para este problema es x
for i = 1:100
for j=2:99
for k = 1:10
evaluaciones(i,j) = evaluaciones(i,j) + (1/(k*pi))*2*((-1)^k)*exp(-cct*(k^2)*(pi^2)*(t(i)))*sin(k*pi*(x(j)));
end
evaluaciones(i,j) = evaluaciones(i,j)+x(j);
end
end
%representación gráfica
surf(t,x,evaluaciones')
title('Solución del problema original')
xlabel('Tiempo')
ylabel('Espacio')
zlabel('Temperatura (°C)')Podemos observar que a medida que aumenta t, la representación de la solución se va pareciendo más a la estacionaria.
A continuación observamos en una gráfica el flujo de calor saliente y entrante en ambos extremos de la varilla a lo largo del tiempo:
YAN
3 Cambio en el coeficiente de conductividad térmica
En este apartado vamos a variar el coeficiente de conductividad térmica de 1 a 1/2. Esto nos deja la siguiente ecuación diferencial:
Nótese que este cambio en el coeficiente de conductividad no tiene influencia en la solución del problema estacionario, por lo que tras homogeneizar las condiciones frontera resulta el siguiente sistema:
Por tanto, la solución para este caso es:
Y su representación gráfica, obtenida con el mismo código que en el apartado anterior pero con coeficiente de conductividad térmica (cct) 1/2 es la siguiente:
Vemos que a simple vista tiene una apariencia casi idéntica a la solución del problema original. Por ello, para entender la influencia del coeficiente de conductividad en la función solución del sistema, veamos las diferencias entre la soluciones [math] u_{1/2} [/math] , [math] u [/math] y la solución estacionaria evaluada en [math] x=\frac{1}{2} [/math]
%Resumen: En este programa, primero creamos vectores de las
%soluciones en x=1/2 con distintos coeficientes de conductividad mediante
% un bucle for. Más tarde se calculan las diferencias en x=1/2 entre estas
%soluciones y la estacionaria, y el %programa acaba gráficando los resultados.
%Definimos la matriz de evaluaciones en x=1/2, y los vectores tiempo y espacio.
%vector de coeficientes de conductividad termica
coefs = [1;1/2;4;1/5;1/10];
evaluaciones= zeros(5,100); %solo las de x=1/2 para cada coeficiente
T = 1;
t = linspace(0,T,100);
x = linspace(0,1,100);
%calculamos las soluciones en x=1/2 para cada coeficiente de coefs
for i = 1:length(coefs)
for k = 1:10
for j = 1:100
evaluaciones(i,j) = evaluaciones(i,j) + (1/(k*pi))*2*((-1)^k)*exp(-coefs(i)*(k^2)*(pi^2)*(t(j)))*sin(k*pi*(x(50)));
end
end
evaluaciones(i,:) = evaluaciones(i,:) + x(50)*ones(1,100); %sumamos la estacionaria
end
%Aquí el vector en x=1/2 de la solución estacionaria.
evaluacioncte=zeros(1,100);
for i = 1:100
evaluacioncte(:,i)=evaluacioncte(:,i)+x(50);
end
%para cada coeficiente graficamos la diferencia
for c = 1:length(coefs)
diferencias = evaluaciones(c,:)-evaluacioncte;
figure(c)
plot(t,diferencias)
title("Diferencia u(1/2,t)-V(1/2) para cct ="+ num2str(coefs(c)) )
xlabel('Tiempo')
ylabel('Diferencia')
ylim([-0.6,0.1])
end
En estas gráficas de diferencias podemos ver como ambas se acercan al valor cero según se va aumentando la [math] t [/math] ya que las soluciones se acercan más a la estacionaria. Además, podemos ver que la solución con el coeficiente de conductividad más bajo lo hace de una forma más lenta que la original. De hecho, probando con distintos valores para este coeficiente podemos concluir que cuanto menor es, más tarda en alcanzarse la solución estacionaria. A continuación se muestran las gráficas obtenidas para coeficientes [math] 4, \frac{1}{5} y \frac{1}{10} [/math] en las que apreciamos mejor el efecto mencionado.
4 Cambio en la condición inicial
Pasamos ahora a variar la condición inicial. Para ello, tendremos que la temperatura en el extremo derecho es también 0 °C, pero consideramos que la temperatura inicial de la barra viene dada por la función [math] u_0(x) = \max\{0, 1-4·|x-1/2|\} [/math]. Así, el sistema EDPS resultante es el siguiente:
En este caso, la solución estacionaria es [math] V(x) = 0 [/math] y por tanto, la gráfica de la temperatura según el espacio y el tiempo es:
%matriz
evaluaciones = zeros(100,100);
%representación gráfica
surf(t,x,evaluaciones')
title('Solución estacionaria')
xlabel('Tiempo')
xlim([0,1])
ylabel('Espacio')
ylim([0,1])
zlabel('Temperatura (°C)')
zlim([-0.1,1])
La solución del sistema por separación de las variables pierde el término x por ser nula la estacionaria y se queda en:
donde [math]c_k = 2·\int_{0}^{1}\sin(k \pi x)·u_0(x)[/math], tras haber considerado la extensión impar de [math]u_0(x)[/math] en [math][-1,1][/math]
Para este sistema no hemos calculado una expresión de la solución, sino que hemos ido calculando los coeficientes de Fourier con la regla del trapecio para cada k y hemos obtenido la siguiente gráfica:
%Resumen: la idea del programa es crear una matriz de pares (tiempo,
%espacio) y en cada uno de esos puntos calcular la aproximación de la
%solución mediante un bucle for que realiza la suma de términos de la serie
%calculada
%creamos la malla en la que introduciremos las evaluaciones
evaluaciones2 = zeros(100,100);
%quedan impuestas los condiciones frontera por ser una matriz de ceros ya
% que posteriormente no modificamos ni la primera ni la última columna
%limite temporal y vectores de tiempos y espacios
T = 1;
t = linspace(0,T,100);
x = linspace(0,1,100);
%funcion de la condicion inicial
cond =@(y) max(0,1-4*abs(y-1/2));
%imponemos la condicion inicial
evaluaciones2(1,:) = cond(x)';
%en el siguiente bucle for, para cada k calculamos el coeficiente de
%Fourier y sumamos el correspondiente término de la serie en cada punto de
% la malla. Para el cálculo de la integral que nos da el coeficiente de
% Fourier usamos la regla del trapecio con 1000 puntos en [0,1]
for k=1:100
%puntos
N=1000;
%extremos del intervalo
a=0; b=1;
h=(b-a)/N;
%partición del intervalo
u=a:h:b;
%f es la funcion que queremos integrar
f=(cond(u).*sin(k.*pi.*u))';
%vector de pesos
w=ones(N+1,1);
w(1)=1/2; w(N+1)=1/2;
%m es el valor de la integral
m=h*w'*f;
%calculo final de los coeficientes del seno
coef = 2*m;
%sumamos los 100 primeros términos de la serie en cada punto de la
% malla con los nuevos coeficientes
for j =2:99
for i = 2:100
evaluaciones2(i,j) = evaluaciones2(i,j) + coef.*sin(k.*pi.*(x(j))).*exp(-(k^2).*(pi^2).*(t(i)));
end
end
end
%plot en 3D
surf(t,x,evaluaciones2')
xlabel('Tiempo')
ylabel('Espacio')
zlabel('Temperatura (°C)')
APARTADO DE LOS FLUJOS - YAN
5 Cambios en las condiciones de frontera y principio del máximo
Por último, repetiremos el apartado anterior pero esta vez suponiendo que la barra está aislada térmicamente en el extremo derecho. Esto quiere decir que pasamos a tener una condicion de frontera de tipo Neumann en [math] x = 1 [/math] :
Si introducimos esta condición en el sistema queda así:
Al resolver por separación de variables obtenemos las siguientes autofunciones:
Véase que esta no es la base trigonométrica usual, y de hecho este conjunto de funciones no tendría por qué ser una base. Sin embargo, se tiene que estas forman un conjunto completo en [0,1] y además se puede comprobar que son ortogonales analíticamente puesto que se cumple que:
Así, podemos concluir que este conjunto de autofunciones es efectivamente una base. Por ello, la solución de nuestro sistema EDPS resulta ser:
donde [math]c_k = 2·\int_{0}^{1}\sin((2k-1)\frac{\pi}{2}x)·u_0(x)[/math], tras haber considerado la extensión impar de [math]u_0(x)[/math] en [math][-1,1][/math]
Gráficamente:
%Resumen: la idea del programa es crear una matriz de pares (tiempo,
%espacio) y en cada uno de esos puntos calcular la aproximación de la
%solución mediante un bucle for que realiza la suma de términos de la serie
%calculada
%creamos la malla en la que introduciremos las evaluaciones
evaluaciones2 = zeros(100,100);
%queda impuesta la condición frontera por ser una matriz de ceros ya
% que posteriormente no modificamos la primera columna
%limite temporal y vectores de tiempos y espacios
T = 1;
t = linspace(0,T,100);
x = linspace(0,1,100);
%funcion de la condicion inicial
cond =@(y) max(0,1-4*abs(y-1/2));
%imponemos la condicion inicial
evaluaciones2(1,:) = cond(x)';
%en el siguiente bucle for, para cada k calculamos el coeficiente de
%Fourier y sumamos el correspondiente término de la serie en cada punto de
% la malla. Para el cálculo de la integral que nos da el coeficiente de
% Fourier usamos la regla del trapecio con 1000 puntos en [0,1]
for k=1:100
%puntos
N=1000;
%extremos del intervalo
a=0; b=1;
h=(b-a)/N;
%partición del intervalo
u=a:h:b;
%f es la funcion que queremos integrar
f=(cond(u).*sin((2*k-1).*(pi/2).*u))';
%vector de pesos
w=ones(N+1,1);
w(1)=1/2; w(N+1)=1/2;
%m es el valor de la integral
m=h*w'*f;
%calculo final de los coeficientes del seno
coef = 2*m;
%sumamos los 100 primeros términos de la serie en cada punto de la
% malla con los nuevos coeficientes
for j =2:100
for i = 2:100
evaluaciones2(i,j) = evaluaciones2(i,j) + coef.*sin((2*k-1).*(pi/2).*(x(j))).*exp(-((2*k-1)^2).*(pi^2)/4.*(t(i)));
end
end
end
%plot en 3D
surf(t,x,evaluaciones2')
xlabel('Tiempo')
ylabel('Espacio')
zlabel('Temperatura (°C)')
Observamos que la nueva gráfica es muy parecida a la anterior y prácticamente lo único que cambia es el extremo [math]x=1[/math], donde hemos impuesto la condición frontera tipo Neumann. En esta gráfica se dibuja una curva en dicho extremo que primero crece un poco más que la del otro extremo y después tiende a la estacionaria (0).
Por último, comprobaremos si se cumple el Principio del Máximo en este problema. Recordamos el enunciado:
[math] \textbf{Teorema: Principio del Máximo.} [/math]
Sea [math] Q_T=\Omega \times (0,T)[/math] con [math] \Omega [/math] abierto de [math] \mathbb{R}^n [/math] y [math] T\gt0 [/math]. Sea [math] u\in C^{2,1}(Q_T)\cap (\bar{Q}_T)[/math] tal que verifica [math] u_t-Δu \leq 0 [/math] en [math] Q_T [/math]. Entonces [math] u [/math] alcanza su máximo en la frontera parabólica [math] \partial _p Q_T [/math], esto es:
6 Ecuación del calor II
En esta sección vamos a dibujar diferentes soluciones de la ecuación del calor usando la solución fundamental.
En primer lugar, tomamos [math]x \in [-1,1] [/math] y [math]t \in [10^{-2},1] [/math] para dibujar la solución fundamental de la ecuación del calor en dimensión 1, que como vimos en clase es:
%Resumen: creamos una matriz que representa una malla de puntos (tiempo,
%espacio) y en cada punto introducimos el valor de la solución fundamental
%evaluada en este
%matriz
evaluaciones = zeros(100,100);
%límite temporal y vectores
T = 1;
t = linspace(10^(-2),T,100);
x = linspace(-1,1,100);
%introducimos los valores de la solucion en cada punto
for i = 1:100
for j=1:100
evaluaciones(i,j) = 1/(2*sqrt(pi*t(i)))*exp(-(x(j)^2)/(4*t(i)));
end
end
%representacion grafica
surf(t,x,evaluaciones')
title('Solución fundamental')
xlabel('Tiempo')
ylabel('Espacio')
zlabel('Temperatura (°C)')
Ahora, vamos a dibujar la solución de la ecuación del calor en una dimensión en el semiespacio [math] x \gt 0 [/math] con condiciones:
Para ello, ...
A contuación veremos la relación de la ecuación del calor con la convolución.
Si tomamos la función [math]u_0(x) = 1_{[-1,1]} [/math] como condición inicial y denotamos la solución fundamental como [math] K(x,t) [/math] , entonces:
es la solución de [math] u_t - u_{xx} = 0[/math] para [math] x \in \mathbb{R} [/math]
Si representamos esta solución para distintos valores de t obtenemos las siguientes gráficas:
Podemos concluir que ...
Para terminar, tomamos el intervalo espacial [math] (x_1,x_2) = [-1,1]^2[/math] para representar la solución fundamental de la ecuación del calor en dimensión 2 para los mismos valores de t:
Observamos que cuanto más cerca está t de 0, más singular es la solución.
