Diferencia entre revisiones de «Ecuación del calor PDM»

De MateWiki
Saltar a: navegación, buscar
Línea 10: Línea 10:
  
  
 +
===Códigos===
 +
============================Decaimiento puntual y en L^2 en caso no acotado========================================
 +
clear; clc; close all;
  
 +
%% Mallado espacial y temporal
 +
x = linspace(-8,8,600);
 +
t = linspace(0.001,10,120);
 +
y = linspace(-1,1,400);  % variable de integracion del dato inicial
  
 +
%% Solucion fundamental
 +
G = @(z,t) (1./sqrt(4*pi*t)) .* exp(-(z.^2)./(4*t));
  
 +
%% Calculo de la solucion U(t,x)
 +
U = zeros(length(t), length(x));
 +
 +
for n = 1:length(t)
 +
    tn = t(n);
 +
    for i = 1:length(x)
 +
        integrando = G(x(i)-y, tn);  % porque u0(y)=1 en [-1,1]
 +
        U(n,i) = trapz(y, integrando);
 +
    end
 +
end
 +
 +
%% =========================
 +
%% 1. Decaimiento puntual
 +
%% =========================
 +
% Elegimos algunos puntos fijos
 +
x_pts = [-0.5, 0, 0.5 1, 2, 4];
 +
idx = zeros(size(x_pts));
 +
 +
for k = 1:length(x_pts)
 +
    [~, idx(k)] = min(abs(x - x_pts(k)));
 +
end
 +
 +
figure;
 +
hold on;
 +
for k = 1:length(x_pts)
 +
    plot(t, U(:,idx(k)), 'LineWidth', 2);
 +
end
 +
grid on;
 +
xlabel('t');
 +
ylabel('u(x_0,t)');
 +
title('Decaimiento puntual de la solucion');
 +
legend('x=-0.5','x=0','x=0.5','x=1','x=2','x=4','Location','best');
 +
 +
%% =========================
 +
%% 2. Decaimiento en norma L2
 +
%% =========================
 +
L2 = zeros(size(t));
 +
 +
for n = 1:length(t)
 +
    L2(n) = sqrt(trapz(x, U(n,:).^2));
 +
end
 +
 +
figure;
 +
plot(t, L2, 'LineWidth', 2);
 +
grid on;
 +
xlabel('t');
 +
ylabel('||u(\cdot,t)||_{L^2}');
 +
title('Decaimiento en norma L^2');
 +
 +
%% =========================
 +
%% 3. Superficie 3D opcional
 +
%% =========================
 +
[X,T] = meshgrid(x,t);
 +
 +
figure;
 +
surf(X,T,U);
 +
shading interp;
 +
colorbar;
 +
xlabel('x');
 +
ylabel('t');
 +
zlabel('u(x,t)');
 +
title('Solucion u(x,t)');
 +
view(135,30);
 +
 +
============================Principio del máximo y aproximación de la solución por convolución========================================
 +
 +
 +
%% PRINCIPIO DEL MAXIMO - ECUACION DEL CALOR EN R
 +
% Casos:
 +
% 1) u0(x) = 1_{[-1,1]}(x)
 +
% 2) u0(x) = exp(-x^2)
 +
%
 +
% Se usa la formula de convolucion:
 +
% u(x,t) = \int Phi(x-y,t) u0(y) dy
 +
% donde Phi(x,t) = 1/sqrt(4*pi*t) * exp(-x^2/(4*t))
 +
 +
clear; close all; clc;
 +
 +
%% Mallado espacial para representar la solucion
 +
x = linspace(-6,6,1600);
 +
 +
%% Instantes de tiempo
 +
tvals = [0.005 0.02 0.1 0.5];
 +
 +
%% Nucleo de calor
 +
Phi = @(z,t) (1./sqrt(4*pi*t)) .* exp(-(z.^2)./(4*t));
 +
 +
%% ============================================================
 +
%% CASO 1: u0(x) = 1_{[-1,1]}(x)
 +
%% ============================================================
 +
 +
u0_1 = @(x) double(abs(x) <= 1);
 +
 +
% Para la integral solo hace falta integrar en [-1,1]
 +
y1 = linspace(-1,1,2500);
 +
 +
figure('Name','Caso 1: u0 = 1_{[-1,1]}');
 +
hold on; grid on;
 +
 +
% Dato inicial
 +
plot(x,u0_1(x),'k--','LineWidth',1.8,'DisplayName','$u_0(x)$');
 +
 +
fprintf('============================\n');
 +
fprintf('CASO 1: u0(x) = 1_{[-1,1]}(x)\n');
 +
fprintf('Dato inicial: min = %.6f, max = %.6f\n', min(u0_1(x)), max(u0_1(x)));
 +
 +
for k = 1:length(tvals)
 +
    t = tvals(k);
 +
 +
    % Matriz Z(i,j)=x_i-y_j
 +
    Z = x(:) - y1(:).';
 +
    integrando = Phi(Z,t);         
 +
    U = trapz(y1, integrando, 2).'; 
 +
 +
    plot(x,U,'LineWidth',1.8,'DisplayName',['$t=', num2str(t), '$']);
 +
 +
    fprintf('t = %.4f --> min = %.6f, max = %.6f\n', t, min(U), max(U));
 +
end
 +
 +
title('Principio del máximo en el primer problema','Interpreter','latex');
 +
xlabel('$x$','Interpreter','latex');
 +
ylabel('$u(x,t)$','Interpreter','latex');
 +
legend('Location','best','Interpreter','latex');
 +
hold off;
 +
 +
%% ============================================================
 +
%% CASO 2: u0(x) = exp(-x^2)
 +
%% ============================================================
 +
 +
u0_2 = @(x) exp(-x.^2);
 +
 +
% Aproximamos la integral en un intervalo grande
 +
y2 = linspace(-8,8,4000);
 +
 +
figure('Name','Caso 2: u0 = exp(-x^2)');
 +
hold on; grid on;
 +
 +
% Dato inicial
 +
plot(x,u0_2(x),'k--','LineWidth',1.8,'DisplayName','$u_0(x)$');
 +
 +
fprintf('\n============================\n');
 +
fprintf('CASO 2: u0(x) = exp(-x^2)\n');
 +
fprintf('Dato inicial: min = %.6f, max = %.6f\n', min(u0_2(x)), max(u0_2(x)));
 +
 +
for k = 1:length(tvals)
 +
    t = tvals(k);
 +
 +
    Z = x(:) - y2(:).';
 +
    integrando = Phi(Z,t) .* exp(-(y2.^2));
 +
    U = trapz(y2, integrando, 2).';
 +
 +
    plot(x,U,'LineWidth',1.8,'DisplayName',['$t=', num2str(t), '$']);
 +
 +
    fprintf('t = %.4f --> min = %.6f, max = %.6f\n', t, min(U), max(U));
 +
end
 +
 +
title('Principio del máximo en el segundo problema','Interpreter','latex');
 +
xlabel('$x$','Interpreter','latex');
 +
 +
ylabel('$u(x,t)$','Interpreter','latex');
 +
legend('Location','best','Interpreter','latex');
 +
hold off;
 +
 +
%% ============================================================
 +
%% EVOLUCION DE MAXIMOS Y MINIMOS
 +
%% ============================================================
 +
 +
tgrid = linspace(0.005,0.5,60);
 +
 +
max1 = zeros(size(tgrid));
 +
min1 = zeros(size(tgrid));
 +
max2 = zeros(size(tgrid));
 +
min2 = zeros(size(tgrid));
 +
 +
for n = 1:length(tgrid)
 +
    t = tgrid(n);
 +
 +
    % Caso 1
 +
    Z1 = x(:) - y1(:).';
 +
    U1 = trapz(y1, Phi(Z1,t), 2).';
 +
    max1(n) = max(U1);
 +
    min1(n) = min(U1);
 +
 +
    % Caso 2
 +
    Z2 = x(:) - y2(:).';
 +
    U2 = trapz(y2, Phi(Z2,t).*exp(-(y2.^2)), 2).';
 +
    max2(n) = max(U2);
 +
    min2(n) = min(U2);
 +
end
 +
 +
figure('Name','Evolucion de maximos y minimos');
 +
 +
subplot(1,2,1)
 +
plot(tgrid,max1,'LineWidth',1.8); hold on; grid on;
 +
plot(tgrid,min1,'LineWidth',1.8);
 +
title('Caso 1: máximos y mínimos','Interpreter','latex');
 +
xlabel('$t$','Interpreter','latex');
 +
ylabel('valor','Interpreter','latex');
 +
legend({'$\max u(\cdot,t)$','$\min u(\cdot,t)$'},'Interpreter','latex','Location','best');
 +
 +
subplot(1,2,2)
 +
plot(tgrid,max2,'LineWidth',1.8); hold on; grid on;
 +
plot(tgrid,min2,'LineWidth',1.8);
 +
title('Caso 2: máximos y mínimos','Interpreter','latex');
 +
xlabel('$t$','Interpreter','latex');
 +
ylabel('valor','Interpreter','latex');
 +
legend({'$\max u(\cdot,t)$','$\min u(\cdot,t)$'},'Interpreter','latex','Location','best');
 +
 +
============================Decaimiento y principio del máximo Dirichlet========================================
 +
clear; clc; close all;
 +
 +
%% ============================================================
 +
%  ECUACION DEL CALOR 1D CON DIRICHLET
 +
%
 +
%    u_t - u_xx = 0,    x in (-6,6), t>0
 +
%    u(-6,t) = 0,        u(6,t) = 0
 +
%    u(x,0) = 1_{[-1,1]}(x)
 +
%
 +
%  Este script:
 +
%  1) Resuelve el problema
 +
%  2) Grafica u(x_i,t) en puntos fijos (decaimiento temporal)
 +
%  3) Muestra la superficie 3D u(x,t) (principio del maximo)
 +
%  4) Grafica perfiles espaciales x -> u(x,t) para tiempos fijos
 +
%% ============================================================
 +
 +
%% Parametros del dominio
 +
a = -6;
 +
b = 6;
 +
Nx = 401;                      % nodos espaciales totales
 +
x = linspace(a,b,Nx)';        % columna
 +
dx = x(2)-x(1);
 +
 +
%% Parametros temporales
 +
Tmax = 8;                      % tiempo final
 +
Nt = 300;                      % numero de tiempos para guardar
 +
tspan = linspace(0,Tmax,Nt);
 +
 +
%% Condicion inicial: 1_{[-1,1]}(x)
 +
u0 = double(abs(x) <= 1);
 +
 +
% Compatibilidad con Dirichlet en bordes
 +
u0(1) = 0;
 +
u0(end) = 0;
 +
 +
u0_max = max(u0);
 +
u0_min = min(u0);
 +
 +
fprintf('Maximo inicial = %.6f\n', u0_max);
 +
fprintf('Minimo inicial = %.6f\n', u0_min);
 +
 +
%% ------------------------------------------------------------
 +
%  Discretizacion espacial
 +
%  Solo resolvemos en nodos interiores
 +
%% ------------------------------------------------------------
 +
Nint = Nx - 2;                % numero de nodos interiores
 +
xint = x(2:end-1);
 +
u0int = u0(2:end-1);
 +
 +
e = ones(Nint,1);
 +
A = spdiags([e -2*e e], -1:1, Nint, Nint) / dx^2;
 +
 +
%% ------------------------------------------------------------
 +
%  Sistema semidiscreto:
 +
%      U_t = A U
 +
%% ------------------------------------------------------------
 +
f = @(t,u) A*u;
 +
 +
% Solver rigido adecuado para difusion
 +
[t,Uint] = ode15s(f, tspan, u0int);
 +
 +
%% Reconstruccion de la solucion completa incluyendo bordes
 +
% Ufull(j,i) = u(x_i, t_j)
 +
Ufull = zeros(length(t), Nx);
 +
Ufull(:,2:end-1) = Uint;
 +
Ufull(:,1) = 0;
 +
Ufull(:,end) = 0;
 +
 +
%% ------------------------------------------------------------
 +
%  Comprobacion numerica del principio del maximo
 +
%% ------------------------------------------------------------
 +
Umax = max(Ufull(:));
 +
Umin = min(Ufull(:));
 +
 +
fprintf('Maximo global numerico de la solucion = %.6f\n', Umax);
 +
fprintf('Minimo global numerico de la solucion = %.6f\n', Umin);
 +
fprintf('Deberia cumplirse aproximadamente: 0 <= u(x,t) <= 1\n');
 +
 +
%% ============================================================
 +
%  1) DECAIMIENTO TEMPORAL EN PUNTOS FIJOS
 +
%% ============================================================
 +
obs_points = [0 0.5 1 2 3 4 5];
 +
obs_idx = zeros(size(obs_points));
 +
obs_real = zeros(size(obs_points));
 +
 +
for k = 1:length(obs_points)
 +
    [~, obs_idx(k)] = min(abs(x - obs_points(k)));
 +
    obs_real(k) = x(obs_idx(k));
 +
end
 +
 +
figure;
 +
hold on;
 +
for k = 1:length(obs_points)
 +
    plot(t, Ufull(:,obs_idx(k)), 'LineWidth', 1.6, ...
 +
        'DisplayName', sprintf('x = %.2f', obs_real(k)));
 +
end
 +
xlabel('t');
 +
ylabel('u(x,t)');
 +
title('Decaimiento temporal en puntos fijos (Dirichlet)');
 +
legend('Location','best');
 +
grid on;
 +
hold off;
 +
 +
 +
%% ============================================================
 +
%  3) PERFILES ESPACIALES x -> u(x,t) PARA TIEMPOS FIJOS
 +
%% ============================================================
 +
times_to_plot = [0 0.19 0.99 3.99 8];
 +
 +
figure;
 +
hold on;
 +
for m = 1:length(times_to_plot)
 +
    [~, j] = min(abs(t - times_to_plot(m)));
 +
    plot(x, Ufull(j,:), 'LineWidth', 1.6, ...
 +
        'DisplayName', sprintf('t = %.2f', t(j)));
 +
end
 +
xlabel('x');
 +
ylabel('u(x,t)');
 +
title('Perfiles espaciales para tiempos fijos');
 +
legend('Location','best');
 +
grid on;
 +
hold off;
 +
 +
============================Decaimiento y principio del máximo Neumann========================================
 +
clear; clc; close all;
 +
 +
%% Dominio espacial
 +
L = 6;
 +
Nx = 401;
 +
x = linspace(-L, L, Nx)';
 +
dx = x(2) - x(1);
 +
 +
%% Tiempo
 +
Tmax = 20;
 +
Nt = 300;
 +
tspan = linspace(0, Tmax, Nt);
 +
 +
%% Condición inicial: 1_{[-1,1]}(x)
 +
u0 = double(abs(x) <= 1);
 +
 +
%% Media integral de la condición inicial
 +
media = trapz(x, u0) / (2*L);
 +
fprintf('Media integral numerica = %.8f\n', media);
 +
fprintf('Media exacta            = %.8f\n', 1/6);
 +
 +
%% Matriz del Laplaciano con Neumann homogéneas
 +
e = ones(Nx,1);
 +
A = spdiags([e -2*e e], -1:1, Nx, Nx);
 +
 +
% Ajuste en bordes para Neumann: u_x = 0
 +
A(1,1)  = -2;
 +
A(1,2)  =  2;
 +
A(end,end-1) = 2;
 +
A(end,end)  = -2;
 +
 +
A = A / dx^2;
 +
 +
%% Resolver sistema semidiscreto U_t = A*U
 +
% ode15s va bien para difusión
 +
f = @(t,u) A*u;
 +
[t,U] = ode15s(f, tspan, u0);
 +
 +
% U sale como matriz Nt x Nx
 +
% para acceder a u(x_i,t_j): U(j,i)
 +
 +
%% Puntos de observación: desde 0 y alejándose
 +
obs_points = [0 0.5 1 2 3 4 5];
 +
obs_idx = zeros(size(obs_points));
 +
obs_real = zeros(size(obs_points));
 +
 +
for k = 1:length(obs_points)
 +
    [~, obs_idx(k)] = min(abs(x - obs_points(k)));
 +
    obs_real(k) = x(obs_idx(k));
 +
end
 +
 +
%% Figura 1: evolución temporal en puntos fijos
 +
figure;
 +
hold on;
 +
for k = 1:length(obs_points)
 +
    plot(t, U(:, obs_idx(k)), 'LineWidth', 1.5, ...
 +
        'DisplayName', sprintf('x = %.2f', obs_real(k)));
 +
end
 +
yline(media, '--k', 'LineWidth', 2, 'DisplayName', 'media = 1/6');
 +
xlabel('t');
 +
ylabel('u(x,t)');
 +
title('Decaimiento temporal en puntos fijos');
 +
legend('Location','best');
 +
grid on;
 +
 +
 +
 +
%% Figura 3: perfiles espaciales para tiempos concretos
 +
times_to_plot = [0 0.2 1 5 20];
 +
figure;
 +
hold on;
 +
for m = 1:length(times_to_plot)
 +
    [~, j] = min(abs(t - times_to_plot(m)));
 +
    plot(x, U(j,:), 'LineWidth', 1.5, ...
 +
        'DisplayName', sprintf('t = %.2f', t(j)));
 +
end
 +
yline(media, '--k', 'LineWidth', 2, 'DisplayName', 'media = 1/6');
 +
xlabel('x');
 +
ylabel('u(x,t)');
 +
title('Perfiles espaciales');
 +
legend('Location','best');
 +
grid on;
  
  
 
[[Categoría:EDP]]
 
[[Categoría:EDP]]
 
[[Categoría:EDP25/26]]
 
[[Categoría:EDP25/26]]

Revisión del 23:32 12 abr 2026

Trabajo realizado por estudiantes
Título Ecuación del calor PDM
Asignatura EDP
Curso 2025-26
Autores Diego García Raposo

Paula Dopico Muñoz

Manuel Herreros Zarco

Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura


1 Ecuación del calor

1.1 Poster

1.2 Códigos

1.2.1 ======================Decaimiento puntual y en L^2 en caso no acotado==================================

clear; clc; close all;

%% Mallado espacial y temporal x = linspace(-8,8,600); t = linspace(0.001,10,120); y = linspace(-1,1,400);  % variable de integracion del dato inicial

%% Solucion fundamental G = @(z,t) (1./sqrt(4*pi*t)) .* exp(-(z.^2)./(4*t));

%% Calculo de la solucion U(t,x) U = zeros(length(t), length(x));

for n = 1:length(t) 

   tn = t(n);
   for i = 1:length(x)
       integrando = G(x(i)-y, tn);   % porque u0(y)=1 en [-1,1]
       U(n,i) = trapz(y, integrando);
   end

end

%% ========================= %% 1. Decaimiento puntual %% ========================= % Elegimos algunos puntos fijos x_pts = [-0.5, 0, 0.5 1, 2, 4]; idx = zeros(size(x_pts));

for k = 1:length(x_pts) 

   [~, idx(k)] = min(abs(x - x_pts(k)));

end

figure; hold on; for k = 1:length(x_pts) 

   plot(t, U(:,idx(k)), 'LineWidth', 2);

end grid on; xlabel('t'); ylabel('u(x_0,t)'); title('Decaimiento puntual de la solucion'); legend('x=-0.5','x=0','x=0.5','x=1','x=2','x=4','Location','best');

%% ========================= %% 2. Decaimiento en norma L2 %% ========================= L2 = zeros(size(t));

for n = 1:length(t) 

   L2(n) = sqrt(trapz(x, U(n,:).^2));

end

figure; plot(t, L2, 'LineWidth', 2); grid on; xlabel('t'); ylabel('||u(\cdot,t)||_{L^2}'); title('Decaimiento en norma L^2');

%% ========================= %% 3. Superficie 3D opcional %% ========================= [X,T] = meshgrid(x,t);

figure; surf(X,T,U); shading interp; colorbar; xlabel('x'); ylabel('t'); zlabel('u(x,t)'); title('Solucion u(x,t)'); view(135,30);

1.2.2 ======================Principio del máximo y aproximación de la solución por convolución==================================

%% PRINCIPIO DEL MAXIMO - ECUACION DEL CALOR EN R % Casos: % 1) u0(x) = 1_{[-1,1]}(x) % 2) u0(x) = exp(-x^2) % % Se usa la formula de convolucion: % u(x,t) = \int Phi(x-y,t) u0(y) dy % donde Phi(x,t) = 1/sqrt(4*pi*t) * exp(-x^2/(4*t))

clear; close all; clc;

%% Mallado espacial para representar la solucion x = linspace(-6,6,1600);

%% Instantes de tiempo tvals = [0.005 0.02 0.1 0.5];

%% Nucleo de calor Phi = @(z,t) (1./sqrt(4*pi*t)) .* exp(-(z.^2)./(4*t));

%% ============================================================ %% CASO 1: u0(x) = 1_{[-1,1]}(x) %% ============================================================

u0_1 = @(x) double(abs(x) <= 1);

% Para la integral solo hace falta integrar en [-1,1] y1 = linspace(-1,1,2500);

figure('Name','Caso 1: u0 = 1_{[-1,1]}'); hold on; grid on;

% Dato inicial plot(x,u0_1(x),'k--','LineWidth',1.8,'DisplayName','$u_0(x)$');

fprintf('============================\n'); fprintf('CASO 1: u0(x) = 1_{[-1,1]}(x)\n'); fprintf('Dato inicial: min = %.6f, max = %.6f\n', min(u0_1(x)), max(u0_1(x)));

for k = 1:length(tvals) 

   t = tvals(k);

   % Matriz Z(i,j)=x_i-y_j
   Z = x(:) - y1(:).';
   integrando = Phi(Z,t);           
   U = trapz(y1, integrando, 2).';  

   plot(x,U,'LineWidth',1.8,'DisplayName',['$t=', num2str(t), '$']);

   fprintf('t = %.4f --> min = %.6f, max = %.6f\n', t, min(U), max(U));

end

title('Principio del máximo en el primer problema','Interpreter','latex'); xlabel('$x$','Interpreter','latex'); ylabel('$u(x,t)$','Interpreter','latex'); legend('Location','best','Interpreter','latex'); hold off;

%% ============================================================ %% CASO 2: u0(x) = exp(-x^2) %% ============================================================

u0_2 = @(x) exp(-x.^2);

% Aproximamos la integral en un intervalo grande y2 = linspace(-8,8,4000);

figure('Name','Caso 2: u0 = exp(-x^2)'); hold on; grid on;

% Dato inicial plot(x,u0_2(x),'k--','LineWidth',1.8,'DisplayName','$u_0(x)$');

fprintf('\n============================\n'); fprintf('CASO 2: u0(x) = exp(-x^2)\n'); fprintf('Dato inicial: min = %.6f, max = %.6f\n', min(u0_2(x)), max(u0_2(x)));

for k = 1:length(tvals) 

   t = tvals(k);

   Z = x(:) - y2(:).';
   integrando = Phi(Z,t) .* exp(-(y2.^2));
   U = trapz(y2, integrando, 2).';

   plot(x,U,'LineWidth',1.8,'DisplayName',['$t=', num2str(t), '$']);

   fprintf('t = %.4f --> min = %.6f, max = %.6f\n', t, min(U), max(U));

end

title('Principio del máximo en el segundo problema','Interpreter','latex'); xlabel('$x$','Interpreter','latex');

ylabel('$u(x,t)$','Interpreter','latex'); legend('Location','best','Interpreter','latex'); hold off;

%% ============================================================ %% EVOLUCION DE MAXIMOS Y MINIMOS %% ============================================================

tgrid = linspace(0.005,0.5,60);

max1 = zeros(size(tgrid)); min1 = zeros(size(tgrid)); max2 = zeros(size(tgrid)); min2 = zeros(size(tgrid));

for n = 1:length(tgrid) 

   t = tgrid(n);

   % Caso 1
   Z1 = x(:) - y1(:).';
   U1 = trapz(y1, Phi(Z1,t), 2).';
   max1(n) = max(U1);
   min1(n) = min(U1);

   % Caso 2
   Z2 = x(:) - y2(:).';
   U2 = trapz(y2, Phi(Z2,t).*exp(-(y2.^2)), 2).';
   max2(n) = max(U2);
   min2(n) = min(U2);

end

figure('Name','Evolucion de maximos y minimos');

subplot(1,2,1) plot(tgrid,max1,'LineWidth',1.8); hold on; grid on; plot(tgrid,min1,'LineWidth',1.8); title('Caso 1: máximos y mínimos','Interpreter','latex'); xlabel('$t$','Interpreter','latex'); ylabel('valor','Interpreter','latex'); legend({'$\max u(\cdot,t)$','$\min u(\cdot,t)$'},'Interpreter','latex','Location','best');

subplot(1,2,2) plot(tgrid,max2,'LineWidth',1.8); hold on; grid on; plot(tgrid,min2,'LineWidth',1.8); title('Caso 2: máximos y mínimos','Interpreter','latex'); xlabel('$t$','Interpreter','latex'); ylabel('valor','Interpreter','latex'); legend({'$\max u(\cdot,t)$','$\min u(\cdot,t)$'},'Interpreter','latex','Location','best');

1.2.3 ======================Decaimiento y principio del máximo Dirichlet==================================

clear; clc; close all;

%% ============================================================ % ECUACION DEL CALOR 1D CON DIRICHLET % % u_t - u_xx = 0, x in (-6,6), t>0 % u(-6,t) = 0, u(6,t) = 0 % u(x,0) = 1_{[-1,1]}(x) % % Este script: % 1) Resuelve el problema % 2) Grafica u(x_i,t) en puntos fijos (decaimiento temporal) % 3) Muestra la superficie 3D u(x,t) (principio del maximo) % 4) Grafica perfiles espaciales x -> u(x,t) para tiempos fijos %% ============================================================

%% Parametros del dominio a = -6; b = 6; Nx = 401;  % nodos espaciales totales x = linspace(a,b,Nx)';  % columna dx = x(2)-x(1);

%% Parametros temporales Tmax = 8;  % tiempo final Nt = 300;  % numero de tiempos para guardar tspan = linspace(0,Tmax,Nt);

%% Condicion inicial: 1_{[-1,1]}(x) u0 = double(abs(x) <= 1);

% Compatibilidad con Dirichlet en bordes u0(1) = 0; u0(end) = 0;

u0_max = max(u0); u0_min = min(u0);

fprintf('Maximo inicial = %.6f\n', u0_max); fprintf('Minimo inicial = %.6f\n', u0_min);

%% ------------------------------------------------------------ % Discretizacion espacial % Solo resolvemos en nodos interiores %% ------------------------------------------------------------ Nint = Nx - 2;  % numero de nodos interiores xint = x(2:end-1); u0int = u0(2:end-1);

e = ones(Nint,1); A = spdiags([e -2*e e], -1:1, Nint, Nint) / dx^2;

%% ------------------------------------------------------------ % Sistema semidiscreto: % U_t = A U %% ------------------------------------------------------------ f = @(t,u) A*u;

% Solver rigido adecuado para difusion [t,Uint] = ode15s(f, tspan, u0int);

%% Reconstruccion de la solucion completa incluyendo bordes % Ufull(j,i) = u(x_i, t_j) Ufull = zeros(length(t), Nx); Ufull(:,2:end-1) = Uint; Ufull(:,1) = 0; Ufull(:,end) = 0;

%% ------------------------------------------------------------ % Comprobacion numerica del principio del maximo %% ------------------------------------------------------------ Umax = max(Ufull(:)); Umin = min(Ufull(:));

fprintf('Maximo global numerico de la solucion = %.6f\n', Umax); fprintf('Minimo global numerico de la solucion = %.6f\n', Umin); fprintf('Deberia cumplirse aproximadamente: 0 <= u(x,t) <= 1\n');

%% ============================================================ % 1) DECAIMIENTO TEMPORAL EN PUNTOS FIJOS %% ============================================================ obs_points = [0 0.5 1 2 3 4 5]; obs_idx = zeros(size(obs_points)); obs_real = zeros(size(obs_points));

for k = 1:length(obs_points) 

   [~, obs_idx(k)] = min(abs(x - obs_points(k)));
   obs_real(k) = x(obs_idx(k));

end

figure; hold on; for k = 1:length(obs_points) 

   plot(t, Ufull(:,obs_idx(k)), 'LineWidth', 1.6, ...
       'DisplayName', sprintf('x = %.2f', obs_real(k)));

end xlabel('t'); ylabel('u(x,t)'); title('Decaimiento temporal en puntos fijos (Dirichlet)'); legend('Location','best'); grid on; hold off;


%% ============================================================ % 3) PERFILES ESPACIALES x -> u(x,t) PARA TIEMPOS FIJOS %% ============================================================ times_to_plot = [0 0.19 0.99 3.99 8];

figure; hold on; for m = 1:length(times_to_plot) 

   [~, j] = min(abs(t - times_to_plot(m)));
   plot(x, Ufull(j,:), 'LineWidth', 1.6, ...
       'DisplayName', sprintf('t = %.2f', t(j)));

end xlabel('x'); ylabel('u(x,t)'); title('Perfiles espaciales para tiempos fijos'); legend('Location','best'); grid on; hold off;

1.2.4 ======================Decaimiento y principio del máximo Neumann==================================

clear; clc; close all;

%% Dominio espacial L = 6; Nx = 401; x = linspace(-L, L, Nx)'; dx = x(2) - x(1);

%% Tiempo Tmax = 20; Nt = 300; tspan = linspace(0, Tmax, Nt);

%% Condición inicial: 1_{[-1,1]}(x) u0 = double(abs(x) <= 1);

%% Media integral de la condición inicial media = trapz(x, u0) / (2*L); fprintf('Media integral numerica = %.8f\n', media); fprintf('Media exacta = %.8f\n', 1/6);

%% Matriz del Laplaciano con Neumann homogéneas e = ones(Nx,1); A = spdiags([e -2*e e], -1:1, Nx, Nx);

% Ajuste en bordes para Neumann: u_x = 0 A(1,1) = -2; A(1,2) = 2; A(end,end-1) = 2; A(end,end) = -2;

A = A / dx^2;

%% Resolver sistema semidiscreto U_t = A*U % ode15s va bien para difusión f = @(t,u) A*u; [t,U] = ode15s(f, tspan, u0);

% U sale como matriz Nt x Nx % para acceder a u(x_i,t_j): U(j,i)

%% Puntos de observación: desde 0 y alejándose obs_points = [0 0.5 1 2 3 4 5]; obs_idx = zeros(size(obs_points)); obs_real = zeros(size(obs_points));

for k = 1:length(obs_points) 

   [~, obs_idx(k)] = min(abs(x - obs_points(k)));
   obs_real(k) = x(obs_idx(k));

end

%% Figura 1: evolución temporal en puntos fijos figure; hold on; for k = 1:length(obs_points) 

   plot(t, U(:, obs_idx(k)), 'LineWidth', 1.5, ...
       'DisplayName', sprintf('x = %.2f', obs_real(k)));

end yline(media, '--k', 'LineWidth', 2, 'DisplayName', 'media = 1/6'); xlabel('t'); ylabel('u(x,t)'); title('Decaimiento temporal en puntos fijos'); legend('Location','best'); grid on;


%% Figura 3: perfiles espaciales para tiempos concretos times_to_plot = [0 0.2 1 5 20]; figure; hold on; for m = 1:length(times_to_plot) 

   [~, j] = min(abs(t - times_to_plot(m)));
   plot(x, U(j,:), 'LineWidth', 1.5, ...
       'DisplayName', sprintf('t = %.2f', t(j)));

end yline(media, '--k', 'LineWidth', 2, 'DisplayName', 'media = 1/6'); xlabel('x'); ylabel('u(x,t)'); title('Perfiles espaciales'); legend('Location','best'); grid on;

Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Ecuación_del_calor_PDM&oldid=104595»