Ver el código de «Ecuación del calor PNL»

{{ TrabajoED | Ecuación del calor. Grupo PNL | [[:Categoría:EDP|EDP]]|[[:Categoría:EDP25/26|2025-26]] | Paula León

Noé Rico

Leo Zambrano  }}

[[Archivo:Poster_PNL2.png||900px]]
[[Archivo:Poster_PNL2.pdf]]

__TOC__

=SOLUCIÓN DEL SISTEMA INICIAL=
El primer código para visualizar la solución del sistema inicial es el siguiente:
<source lang: "Matlab" line>

% Parametros
x = linspace(0, 1, 1000);    % Mallado espacial
t = linspace(0, 1, 1000);    % Mallado temporal
[X, T] = meshgrid(x, t);    % Malla 2D

% Solución estacionaria
v = 9*X + 1;

% Primeros terminos de Fourier
N = 10;                     % Número de términos
w = zeros(size(X));

for n = 1:N
    cn = 18 / (n * pi);
    wn = cn * sin(n * pi .* X) .* exp(-(n * pi)^2 .* T);
    w = w + wn;
end

% Solución tras el cambio
u = w + v;

% Gráfica
figure;
surf(X, T, u, 'EdgeColor', 'none');
xlabel('x');
ylabel('t');
zlabel('u(x,t)');
title('Solución de la ecuación del calor: u(x,t) con primeros 10 términos');
colormap('jet');
colorbar;
view(45, 30);   
grid on;

</source>

=FLUJO EN LOS EXTREMOS DEL SI=
El código para observar el flujo en los extremos es:

<source lang: "Matlab" line>

% Parámetros
t = linspace(1e-16, 1, 1000);   % tiempo, evitando t=0 exacto
N = 10;                         % número de términos

% Inicializar flujos
Phi0 = -9 * ones(size(t));
Phi1 = -9 * ones(size(t));

% Sumar los N primeros términos
for n = 1:N
    Phi0 = Phi0 - 18 * exp(-(n*pi)^2 .* t);
    Phi1 = Phi1 - 18 * (-1)^n * exp(-(n*pi)^2 .* t);
end

% Gráfica
figure;
plot(t, Phi0, 'b-', 'LineWidth', 1.5); hold on;
plot(t, Phi1, 'r-', 'LineWidth', 1.5);
xlabel('tiempo t');
ylabel('Flujo de calor Φ(x,t)');
legend('Φ(0,t) (extremo izquierdo)', 'Φ(1,t) (extremo derecho)', 'Location', 'best');
title('Flujo de calor en los extremos de la varilla');
grid on;

</source>

=CAMBIO DE COEFICIENTE DE DIFUSIÓN=
El código utilizado para la solución con coeficiente de difusión 0.5 y la comparación con el coeficiente de difusión 1 es el siguiente:

<source lang: "Matlab" line>

% Parámetros
x = linspace(0, 1, 1000);      % mallado espacial
t = linspace(0, 1, 1000);      % mallado temporal (incluye t=0)
[X, T] = meshgrid(x, t);      % malla 2D

% Solución estacionaria 
v = 9*X + 1;                  

% Número de términos de la serie
N = 10;

% Coeficientes c_n 
c = zeros(1, N);
for n = 1:N
    c(n) = 18 / (n * pi);
end

% 1. Solución para k = 0.5
k1 = 0.5;
w_k05 = zeros(size(X));
for n = 1:N
    wn = c(n) * sin(n*pi*X) .* exp(-k1 * (n*pi)^2 .* T);
    w_k05 = w_k05 + wn;
end
u_k05 = v + w_k05;

% Gráfica 3D para k=0.5
figure;
surf(X, T, u_k05, 'EdgeColor', 'none');
xlabel('x'); ylabel('t'); zlabel('u(x,t)');
title('Solución con conductividad k = 0.5 (primeros 10 términos)');
colormap('jet'); colorbar; grid on;


% 2. Diferencia u(1/2,t)-uest(1/2) para ambos k
% Punto medio x=0.5
x_mid = 0.5;
% Para k=1 (original)
k_orig = 1;
w_diff_orig = zeros(size(t));
for n = 1:N
    wn = c(n) * sin(n*pi*x_mid) .* exp(-k_orig * (n*pi)^2 .* t);
    w_diff_orig = w_diff_orig + wn;
end
% Para k=0.5
w_diff_k05 = zeros(size(t));
for n = 1:N
    wn = c(n) * sin(n*pi*x_mid) .* exp(-k1 * (n*pi)^2 .* t);
    w_diff_k05 = w_diff_k05 + wn;
end

% Valor estacionario en x=0.5
uest_mid = 9*x_mid + 1;   % = 5.5
% La diferencia es simplemente w(x=0.5,t)
diff_orig = w_diff_orig;
diff_k05 = w_diff_k05;

% Gráfica comparativa
figure;
plot(t, diff_orig, 'b-', 'LineWidth', 1.5); hold on;
plot(t, diff_k05, 'r--', 'LineWidth', 1.5);
xlabel('tiempo t'); ylabel('u(0.5,t) - u_{est}(0.5)');
legend('k = 1', 'k = 0.5', 'Location', 'best');
title('Evolución de la diferencia respecto al estado estacionario en x = 0.5');
grid on;

</source>

=SOLUCIÓN CON CAMBIO EN LAS CONDICIONES INICIALES=
El código para visualizar la solución del sistema con las condiciones iniciales modificadas es:

<source lang: "Matlab" line>

% Parametros
x = linspace(0, 1, 1000);    % Mallado espacial
t = linspace(0, 1, 1000);    % Mallado temporal
[X, T] = meshgrid(x, t);    % Malla 2D

% Solución estacionaria
v = 10;

% Primeros terminos de Fourier
N = 10;                     % Número de términos
w = zeros(size(X));

for n = 1:N 
    cn = 20/(n*pi)*(cos(2*pi*n/3)-cos(pi*n/3));
    wn = cn * sin(n * pi .* X) .* exp(-(n * pi)^2 .* T);
    w = w + wn;
end

% Solución tras el cambio
u = w + v;

% Gráfica
figure;
surf(X, T, u, 'EdgeColor', 'none');
xlabel('x');
ylabel('t');
zlabel('u(x,t)');
title('Solución de la ecuación del calor: u(x,t) con primeros 10 términos');
colormap('jet');
colorbar;
view(45, 30);   
grid on;

</source>

=FLUJO EN LOS EXTREMOS CON CAMBIO EN CONDICIONES INICIALES=
Para visualizar el flujo con las nuevas condiciones iniciales utilizamos el siguiente código:

<source lang: "Matlab" line>

% Parámetros
t = linspace(1e-16, 1, 1000);   % tiempo, evitando t=0 exacto
N = 10;                         % número de términos

% Inicializar flujos
Phi0 = zeros(size(t));
Phi1 = zeros(size(t));

% Sumar los N primeros términos
for n = 1:N
    Phi0 = Phi0 - 20 * (cos(2*pi*n/3)-cos(pi*n/3))* exp(-(n*pi)^2 .* t);
    Phi1 = Phi1 - 20 * (cos(2*pi*n/3)-cos(pi*n/3))*(-1)^n * exp(-(n*pi)^2 .* t);
end

% Gráfica
figure;
plot(t, Phi0, 'b-', 'LineWidth', 1.5); hold on;
plot(t, Phi1, 'r-', 'LineWidth', 1.5);
xlabel('tiempo t');
ylabel('Flujo de calor Φ(x,t)');
legend('Φ(0,t) (extremo izquierdo)', 'Φ(1,t) (extremo derecho)', 'Location', 'best');
title('Flujo de calor en los extremos de la varilla');
grid on;

</source>

=SOLUCIÓN CON CAMBIO EN LAS CONDICIONES FRONTERA=
El código para observar la solución con el sistema con distintas condiciones de frontera es:

<source lang: "Matlab" line>

% Parametros
x = linspace(0, 1, 1000);    % Mallado espacial
t = linspace(0, 1, 1000);    % Mallado temporal
[X, T] = meshgrid(x, t);    % Malla 2D

% Solución estacionaria
v = 10;

% Primeros terminos de Fourier
N = 10;                     % Número de términos
w = zeros(size(X));

for n = 0:N 
    cn = -40/((n+1/2)*pi)*sin((pi/12) + (n*pi/6)) * cos((pi/4) + (n*pi/2));
    wn = cn * cos((1/2 + n)*pi.*X).*exp(-((1/2+n)*pi)^2.*T);
    w = w + wn;
end

% Solución tras el cambio
u = w + v;

% Gráfica
figure;
surf(X, T, u, 'EdgeColor', 'none');
xlabel('x');
ylabel('t');
zlabel('u(x,t)');
title('Solución de la ecuación del calor: u(x,t) con primeros 10 términos');
colormap('jet');
colorbar;
view(45, 30);   
grid on;

</source>

=FLUJO EN LOS EXTREMOS CON CAMBIO EN CONDICIONES FRONTERA=
Para ver el flujo en los extremos del problema con distintas condiciones de fronteras utilizamos el siguiente código:

<source lang: "Matlab" line>

% Parámetros
t = linspace(1e-16, 1, 1000);   % tiempo, evitando t=0 exacto
N = 10;                         % número de términos

% Inicializar flujos
Phi0 = zeros(size(t));
Phi1 = zeros(size(t));

% Sumar los N primeros términos
for n = 1:N
    Phi0 = Phi0 + 40*sin((pi/12) + (n*pi/6)) * cos((pi/4) + (n*pi/2)) * sin((1/2 + n)*pi.*0)*exp(-((1/2+n)*pi)^2.*t);
    Phi1 = Phi1 + 40*sin((pi/12) + (n*pi/6)) * cos((pi/4) + (n*pi/2)) * sin((1/2 + n)*pi)*exp(-((1/2+n)*pi)^2.*t);
end

% Gráfica
figure;
plot(t, Phi0, 'b-', 'LineWidth', 1.5); hold on;
plot(t, Phi1, 'r-', 'LineWidth', 1.5);
xlabel('tiempo t');
ylabel('Flujo de calor Φ(x,t)');
legend('Φ(0,t) (extremo izquierdo)', 'Φ(1,t) (extremo derecho)', 'Location', 'best');
title('Flujo de calor en los extremos de la varilla');
grid on;

</source>

=PRINCIPIO DE COMPARACIÓN=
El código para visualizar el principio de comparación es:

<source lang: "Matlab" line>

% Parámetros
x = linspace(0, 1, 1000);      % Mallado espacial
t = linspace(0, 1, 1000);      % Mallado temporal
[X, T] = meshgrid(x, t);       % Malla 2D

% Solución estacionaria (cota inferior)
v = 9 * X + 1;                  % v(x,t) = 9x+1

% Cota superior constante
z = 10 * ones(size(X));         % z(x,t) = 10

% Cálculo de la solución u(x,t) del problema original
% u_t = u_xx, con u(0,t)=1, u(1,t)=10, u(x,0)=10
% La solución se escribe como u = v + w, donde w se obtiene por separación de variables
N = 1000;                         % Número de términos de Fourier
w = zeros(size(X));

for n = 1:N
    cn = 18 / (n * pi);         % Coeficientes de la serie de Fourier de 10 - (9x+1) = 9(1-x)
    wn = cn * sin(n * pi .* X) .* exp(-(n * pi)^2 .* T);
    w = w + wn;
end

u = v + w;                      % Solución completa

% Gráfica 3D que muestra u(x,t) y las dos cotas
figure;
% Superficie de la solución
surf(X, T, u, 'EdgeColor', 'none', 'FaceAlpha', 0.8);
hold on;
% Cota inferior v(x,t) = 9x+1
surf(X, T, v, 'EdgeColor', 'none', 'FaceAlpha', 0.4, 'FaceColor', 'green');
% Cota superior z(x,t) = 10
surf(X, T, z, 'EdgeColor', 'none', 'FaceAlpha', 0.4, 'FaceColor', 'red');

xlabel('x');
ylabel('t');
zlabel('u(x,t)');
title('Comparación: solución u(x,t) acotada entre 9x+1 (verde) y 10 (rojo)');
legend('u(x,t)', '9x+1 (cota inferior)', '10 (cota superior)', 'Location', 'best');
colormap('jet');
colorbar;
view(45, 30);   
grid on;
hold off;

</source>

[[Categoría:EDP]]
[[Categoría:EDP25/26]]