Ecuación del calor PDM

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.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.1.1 ======================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.1.2 ======================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.1.3 ======================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;