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

Revisión actual del 08:44 15 abr 2026

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;