Diferencia entre revisiones de «Parte de Andrews y Lucía»

Revisión actual del 20:28 7 mar 2024

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]. CREO QUE AITANA TIENE QUE DECIR COSAS AQUI

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,

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 más marcada que la del flujo saliente

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,

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, ESCRIBIR LA SOLUCIÓN FINAL 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. Ambas gráficas se han dibujado con el código,

INTRODUCIR CÓDIGO

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]. INTRODUCIR CÓDIGO

centro

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,

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,

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]. NO SE QUE DECIRLE AQUI

3.2 Flujo del calor entrante y saliente

Flujos entrante y saliente

NI IDEA

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:

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:

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.1 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.