MateWiki - Contribuciones del usuario [es]

Ecuación del Calor Estacionaria (Grupo DIF)

2025-11-12T10:09:32Z

Daniel Rodríguez Calderón:

{{ TrabajoED | Ecuación del Calor Estacionaria (Grupo DIF) | [[:Categoría:MNEDP|MNEDP]]|[[:Categoría:MNEDP25/26|2025-26]] |
*Daniel Rodríguez Calderón,
*Israel López Moreno,
*Francisco Lavao Amengual.}}

=Póster del Trabajo=

{|style="margin: 0 auto;"
| [[Archivo:Poster_MNEDP_DIF.jpg|900px|thumb|right|Trabajo de Métodos Numéricos en Ecuaciones Diferenciales en Derivadas Parciales sobre la Ecuación del Calor Estacionaria en formato póster realizado por el grupo DIF.]]
|}

=Códigos de Matlab=
<source lang="matlab">
%% MEF de la Ecuación del Calor con condiciones Dirichlet homogéneas + Neumann homogéneas
clc; clear all; close all

%% Ejemplo de Dominio Cuadrado Unidad

%%% --- Parámetros de la malla ---
n_x = 30; n_y = 30; tol = 1e-12;

[x_v, y_v] = meshgrid(linspace(0,1,n_x+1), linspace(0,1,n_y+1));
X = x_v(:); Y = y_v(:);
nodos = [X Y];

triang = delaunay(nodos);
n_nodos = size(nodos,1);
n_triang = size(triang,1);

figure(1)
subplot(1,3,1)
trimesh(triang, nodos(:,1), nodos(:,2), zeros(size(nodos,1),1))
hold on; xlabel('x'), ylabel('y');
%shading interp; camproj perspective; axis vis3d; view(45,30)
title('\textbf{Dominio cuadrado unidad}', 'Interpreter','latex')
scatter(0.1, 0.5, 'red', '*', LineWidth = 3)
view(2); axis equal tight;

%%% --- Parámetros del problema ---
k_fun = @(x,y) 1; % conductividad
c_fun = @(x,y) 0; % término de reacción
% r_fun = @(x,y) 10*exp(-50*((x-0.5).^2+(y-0.5).^2)); % fuente (como en tus códigos)
r_fun = @(x,y) 350*exp(-10*((x-0.1).^2+(y-0.5).^2)); % Fuente Gaussiana

%%% --- Selección de frontera Neumann ---
% Opciones: 'gamma1','gamma2','gamma3','gamma4'
neumann_sides = ["gamma1","gamma2","gamma3"];
g_val = 0; % flujo de Neumann en la pared seleccionada (0 según tu petición)

%%% --- Ensamblaje global ---
I = []; J = []; V = []; R = zeros(n_nodos,1);

for t = 1:n_triang
id = triang(t,:);
puntos = nodos(id,:);
x1 = puntos(1,1); y1 = puntos(1,2);
x2 = puntos(2,1); y2 = puntos(2,2);
x3 = puntos(3,1); y3 = puntos(3,2);

Area = abs(det([x2-x1 x3-x1; y2-y1 y3-y1]))/2;
B = [y2-y3 y3-y1 y1-y2; x3-x2 x1-x3 x2-x1]/(2*Area);
x_m = mean(puntos(:,1)); y_m = mean(puntos(:,2));

k = k_fun(x_m, y_m);
Ke = k*Area*(B'*B);

c = c_fun(x_m, y_m);
Me = (c*Area/12)*[2 1 1; 1 2 1; 1 1 2];

rmid = r_fun(x_m, y_m);
Re = rmid*Area/3*ones(3,1);

for a = 1:3
for b = 1:3
I(end+1) = id(a);
J(end+1) = id(b);
V(end+1) = Ke(a,b) + Me(a,b);
end
R(id(a)) = R(id(a)) + Re(a);
end
end

A = sparse(I, J, V, n_nodos, n_nodos);
b = R;

%%% --- Identificación de nodos de contorno y clasificación Dirichlet/Neumann ---
is_gamma1 = abs(nodos(:,2)) < tol & nodos(:,1) >= 0 & nodos(:,1) <= 1;
is_gamma2 = abs(nodos(:,1)) < tol & nodos(:,2) >= 0 & nodos(:,2) <= 1;
is_gamma3 = abs(nodos(:,2)-1) < tol & nodos(:,1) >= 0 & nodos(:,1) <= 1;
is_gamma4 = abs(nodos(:,1)-1) < tol & nodos(:,2) >= 0 & nodos(:,2) <= 1;

% Por defecto: Dirichlet en todas las paredes
dirichlet_nodes = is_gamma1 | is_gamma2 | is_gamma3 | is_gamma4;
neumann_nodes = zeros(length(is_gamma1), 1);

% Reemplazar la pared seleccionada por Neumann
for i = 1:numel(neumann_sides)
neumann_side = neumann_sides(i);

switch lower(neumann_side)
case 'gamma1'
dirichlet_nodes(is_gamma1) = false;
neumann_nodes = neumann_nodes + is_gamma1;
case 'gamma2'
dirichlet_nodes(is_gamma2) = false;
neumann_nodes = neumann_nodes + is_gamma2;
case 'gamma3'
dirichlet_nodes(is_gamma3) = false;
neumann_nodes = neumann_nodes + is_gamma3;
case 'gamma4'
dirichlet_nodes(is_gamma4) = false;
neumann_nodes = neumann_nodes + is_gamma4;
otherwise
error('neumann_side no reconocido. Usa gamma1,gamma2,gamma3,gamma4,gamma5,gamma6.')
end
end

ind_dir = find(dirichlet_nodes);
ind_neu = find(neumann_nodes);

%%% --- Aplicación de condiciones Dirichlet homogéneas (en nodos dirichlet_nodes) ---
% Igual que en tu primer código: anular filas/columnas y poner 1 en diagonal
A(ind_dir, :) = 0; A(:, ind_dir) = 0;
A(sub2ind(size(A), ind_dir, ind_dir)) = 1;
b(ind_dir) = 0;

%%% --- Resolución del sistema ---
u = A\b;
u = u(:);

%%% --- Visualización ---
figure(2);
subplot(3,3,1)
trisurf(triang, nodos(:,1), nodos(:,2), u, 'EdgeColor', 'none');
view(2); axis equal tight; colorbar; colormap(turbo);
title('\textbf{Solucion FEM (vista superior) sobre el dominio cuadrado unidad}', 'Interpreter','latex')
xlabel('x'), ylabel('y')

subplot(3,3,2)
trisurf(triang, nodos(:,1), nodos(:,2), u)
shading interp; camproj perspective; axis vis3d; view(45,30)
xlabel('x'), ylabel('y'), zlabel('u')
title('\textbf{Solucion FEM (superficie) sobre el dominio cuadrado unidad}', 'Interpreter','latex')

% Mostrar en la malla qué nodos son Dirichlet/Neumann
subplot(3,3,3)
plot(nodos(:,1), nodos(:,2), '.k'); hold on
plot(nodos(ind_dir,1), nodos(ind_dir,2), 'or', 'MarkerFaceColor','r');
plot(nodos(ind_neu,1), nodos(ind_neu,2), 'sb', 'MarkerFaceColor','b');
legend('Interior/otros','Dirichlet','Neumann','Location','bestoutside')
view(2); axis equal tight; xlabel('x'), ylabel('y');
title('\textbf{Nodos de contorno sobre el dominio cuadrado unidad}', 'Interpreter','latex');

%%% --- Barrido de parámetros k (igual que en tus scripts originales) ---
k_vals = [0.1 1 5 10 15 20];
maxU = zeros(size(k_vals));

for i = 1:length(k_vals)
k_fun = @(x,y) k_vals(i);

% Re-ensamblar para cada k
I = []; J = []; V = []; R = zeros(n_nodos,1);
for t = 1:n_triang
id = triang(t,:);
puntos = nodos(id,:);
x1 = puntos(1,1); y1 = puntos(1,2);
x2 = puntos(2,1); y2 = puntos(2,2);
x3 = puntos(3,1); y3 = puntos(3,2);

Area = abs(det([x2-x1 x3-x1; y2-y1 y3-y1]))/2;
B = [y2-y3 y3-y1 y1-y2; x3-x2 x1-x3 x2-x1]/(2*Area);
x_m = mean(puntos(:,1)); y_m = mean(puntos(:,2));

k = k_fun(x_m, y_m);
Ke = k*Area*(B'*B);
c = c_fun(x_m, y_m);
Me = (c*Area/12)*[2 1 1; 1 2 1; 1 1 2];
rmid = r_fun(x_m, y_m);
Re = rmid*Area/3*ones(3,1);

for a = 1:3
for b = 1:3
I(end+1) = id(a);
J(end+1) = id(b);
V(end+1) = Ke(a,b) + Me(a,b);
end
R(id(a)) = R(id(a)) + Re(a);
end
end
A_k = sparse(I, J, V, n_nodos, n_nodos);
b_k = R;

% Impón Dirichlet (mismo conjunto de nodos)
A_k(ind_dir, :) = 0; A_k(:, ind_dir) = 0;
A_k(sub2ind(size(A_k), ind_dir, ind_dir)) = 1;
b_k(ind_dir) = 0;

% Resolver y almacenar max
u_k = A_k \ b_k;
maxU(i) = max(u_k);
end

figure(3)
subplot(1,3,1)
plot(k_vals, maxU, '-o', 'LineWidth', 1.5);
xlabel('$k$', 'Interpreter','latex'); ylabel('\textbf{max \{u\}}', 'Interpreter','latex');
title('\textbf{Efecto de $k$ en la solucion sobre el dominio cuadrado unidad}', 'Interpreter','latex');
grid on; hold on; view(2); axis square;

%% Ejemplo de Dominio Tipo L

%%% --- Parámetros de la malla ---
n_x = 30; n_y = 30;
[x_v, y_v] = meshgrid(linspace(0,3,n_x+1), linspace(0,2,n_y+1));
X = x_v(:); Y = y_v(:);
nodos = [X Y];

fuera_dominio = nodos(:,1) > 1 & nodos(:,2) < 1;
nodos = nodos(~fuera_dominio,:);

tol = 1e-12;

is_gamma1 = abs(nodos(:,2)) < tol & nodos(:,1) >= 0 & nodos(:,1) <= 1;
is_gamma2 = abs(nodos(:,1)) < tol & nodos(:,2) >= 0 & nodos(:,2) <= 2;
is_gamma3 = abs(nodos(:,2)-2) < tol & nodos(:,1) >= 0 & nodos(:,1) <= 3;
is_gamma4 = abs(nodos(:,1)- 3) < tol & nodos(:,2) >= 1 & nodos(:,2) <= 2;
is_gamma5 = abs(nodos(:,2)-1) < tol & nodos(:,1) >= 1 & nodos(:,1) <= 3;
is_gamma6 = abs(nodos(:,1)-1) < tol & nodos(:,2) >= 0 & nodos(:,2) <= 1;
is_beta = nodos(:,1) > 1 & nodos(:,1) <= 3 & nodos(:,2) >= 0 & nodos(:,2) < 1;

triang = delaunay(nodos);
n_nodos = size(nodos,1);

is_edge = is_gamma5 | is_gamma6;
triang_eliminar = all(is_edge(triang), 2);
triang(triang_eliminar, :) = [];

n_triang = size(triang,1);

figure(1)
subplot(1,3,2)
trimesh(triang, nodos(:,1), nodos(:,2), zeros(size(nodos,1),1))
%shading interp; camproj perspective; axis vis3d; view(45,30)
xlabel('x'), ylabel('y'); hold on;
title('\textbf{Dominio tipo $L$}', 'Interpreter','latex')
scatter([0.5 1.5], [0.1 1.9], 'red', '*', LineWidth = 3)
view(2); axis equal tight;

%%% --- Parámetros del problema ---
k_fun = @(x,y) 1; % conductividad
c_fun = @(x,y) 0; % término de reacción
% r_fun = @(x,y) 10*exp(-50*((x-0.5).^2+(y-0.5).^2)); % fuente (como en tus códigos)
r_fun = @(x,y) 75*exp(-10*((x-0.5).^2+(y-0.1).^2)) + 75*exp(-10*((x-1.5).^2+(y-1.9).^2)); % Fuente Gaussiana

%%% --- Selección de frontera Neumann ---
% Opciones: 'gamma1','gamma2','gamma3','gamma4','gamma5','gamma6'
neumann_sides = ["gamma1","gamma2","gamma3","gamma5","gamma6"];
g_val = 0; % flujo de Neumann en la pared seleccionada (0 según tu petición)

%%% --- Ensamblaje global ---
I = []; J = []; V = []; R = zeros(n_nodos,1);

for t = 1:n_triang
id = triang(t,:);
puntos = nodos(id,:);
x1 = puntos(1,1); y1 = puntos(1,2);
x2 = puntos(2,1); y2 = puntos(2,2);
x3 = puntos(3,1); y3 = puntos(3,2);

Area = abs(det([x2-x1 x3-x1; y2-y1 y3-y1]))/2;
B = [y2-y3 y3-y1 y1-y2; x3-x2 x1-x3 x2-x1]/(2*Area);
x_m = mean(puntos(:,1)); y_m = mean(puntos(:,2));

k = k_fun(x_m, y_m);
Ke = k*Area*(B'*B);

c = c_fun(x_m, y_m);
Me = (c*Area/12)*[2 1 1; 1 2 1; 1 1 2];

rmid = r_fun(x_m, y_m);
Re = rmid*Area/3*ones(3,1);

for a = 1:3
for b = 1:3
I(end+1) = id(a);
J(end+1) = id(b);
V(end+1) = Ke(a,b) + Me(a,b);
end
R(id(a)) = R(id(a)) + Re(a);
end
end

A = sparse(I, J, V, n_nodos, n_nodos);
b = R;

%%% --- Identificación de nodos de contorno y clasificación Dirichlet/Neumann ---
is_gamma1 = abs(nodos(:,2)) < tol & nodos(:,1) >= 0 & nodos(:,1) <= 1;
is_gamma2 = abs(nodos(:,1)) < tol & nodos(:,2) >= 0 & nodos(:,2) <= 2;
is_gamma3 = abs(nodos(:,2)-2) < tol & nodos(:,1) >= 0 & nodos(:,1) <= 3;
is_gamma4 = abs(nodos(:,1)- 3) < tol & nodos(:,2) >= 1 & nodos(:,2) <= 2;
is_gamma5 = abs(nodos(:,2)-1) < tol & nodos(:,1) >= 1 & nodos(:,1) <= 3;
is_gamma6 = abs(nodos(:,1)-1) < tol & nodos(:,2) >= 0 & nodos(:,2) <= 1;
is_beta = nodos(:,1) > 1 & nodos(:,1) <= 3 & nodos(:,2) >= 0 & nodos(:,2) < 1;

% Por defecto: Dirichlet en todas las paredes
dirichlet_nodes = is_gamma1 | is_gamma2 | is_gamma3 |...
is_gamma4 | is_gamma5 | is_gamma6 | is_beta;
neumann_nodes = zeros(length(is_gamma1),1);

% Reemplazar la pared seleccionada por Neumann
for i = 1:numel(neumann_sides)
neumann_side = neumann_sides(i);

switch lower(neumann_side)
case 'gamma1'
dirichlet_nodes(is_gamma1) = false;
neumann_nodes = neumann_nodes + is_gamma1;
case 'gamma2'
dirichlet_nodes(is_gamma2) = false;
neumann_nodes = neumann_nodes + is_gamma2;
case 'gamma3'
dirichlet_nodes(is_gamma3) = false;
neumann_nodes = neumann_nodes + is_gamma3;
case 'gamma4'
dirichlet_nodes(is_gamma4) = false;
neumann_nodes = neumann_nodes + is_gamma4;
case 'gamma5'
dirichlet_nodes(is_gamma5) = false;
neumann_nodes = neumann_nodes + is_gamma5;
case 'gamma6'
dirichlet_nodes(is_gamma6) = false;
neumann_nodes = neumann_nodes + is_gamma6;
otherwise
error('neumann_side no reconocido. Usa gamma1,gamma2,gamma3,gamma4,gamma5,gamma6.')
end
end

ind_dir = find(dirichlet_nodes);
ind_neu = find(neumann_nodes);

%%% --- Aplicación de condiciones Dirichlet homogéneas (en nodos dirichlet_nodes) ---
% Igual que en tu primer código: anular filas/columnas y poner 1 en diagonal
A(ind_dir, :) = 0; A(:, ind_dir) = 0;
A(sub2ind(size(A), ind_dir, ind_dir)) = 1;
b(ind_dir) = 0;

%%% --- Resolución del sistema ---
u = A\b;
u = u(:);

%%% --- Visualización ---
figure(2);
subplot(3,3,4)
trisurf(triang, nodos(:,1), nodos(:,2), u, 'EdgeColor', 'none');
view(2); axis equal tight; colorbar;
title('\textbf{Solucion FEM (vista superior) sobre el dominio tipo $L$}', 'Interpreter','latex')
xlabel('x'), ylabel('y')

subplot(3,3,5)
trisurf(triang, nodos(:,1), nodos(:,2), u)
shading interp; camproj perspective; axis vis3d; view(45,30)
xlabel('x'), ylabel('y'), zlabel('u')
title('\textbf{Solucion FEM (superficie) sobre el dominio tipo $L$}', 'Interpreter','latex')

% Mostrar en la malla qué nodos son Dirichlet/Neumann
subplot(3,3,6)
plot(nodos(:,1), nodos(:,2), '.k'); hold on
plot(nodos(ind_dir,1), nodos(ind_dir,2), 'or', 'MarkerFaceColor','r');
plot(nodos(ind_neu,1), nodos(ind_neu,2), 'sb', 'MarkerFaceColor','b');
legend('Interior/otros','Dirichlet','Neumann','Location','bestoutside')
view(2); axis equal tight; xlabel('x'), ylabel('y')
title('\textbf{Nodos de contorno sobre el dominio tipo $L$}', 'Interpreter','latex');

%%% --- Barrido de parámetros k (igual que en tus scripts originales) ---
k_vals = [0.1 1 5 10 15 20];
maxU = zeros(size(k_vals));

for i = 1:length(k_vals)
k_fun = @(x,y) k_vals(i);

% Re-ensamblar para cada k
I = []; J = []; V = []; R = zeros(n_nodos,1);
for t = 1:n_triang
id = triang(t,:);
puntos = nodos(id,:);
x1 = puntos(1,1); y1 = puntos(1,2);
x2 = puntos(2,1); y2 = puntos(2,2);
x3 = puntos(3,1); y3 = puntos(3,2);

Area = abs(det([x2-x1 x3-x1; y2-y1 y3-y1]))/2;
B = [y2-y3 y3-y1 y1-y2; x3-x2 x1-x3 x2-x1]/(2*Area);
x_m = mean(puntos(:,1)); y_m = mean(puntos(:,2));

k = k_fun(x_m, y_m);
Ke = k*Area*(B'*B);
c = c_fun(x_m, y_m);
Me = (c*Area/12)*[2 1 1; 1 2 1; 1 1 2];
rmid = r_fun(x_m, y_m);
Re = rmid*Area/3*ones(3,1);

for a = 1:3
for b = 1:3
I(end+1) = id(a);
J(end+1) = id(b);
V(end+1) = Ke(a,b) + Me(a,b);
end
R(id(a)) = R(id(a)) + Re(a);
end
end
A_k = sparse(I, J, V, n_nodos, n_nodos);
b_k = R;

% Impón Dirichlet (mismo conjunto de nodos)
A_k(ind_dir, :) = 0; A_k(:, ind_dir) = 0;
A_k(sub2ind(size(A_k), ind_dir, ind_dir)) = 1;
b_k(ind_dir) = 0;

% Resolver y almacenar max
u_k = A_k \ b_k;
maxU(i) = max(u_k);
end

figure(3)
subplot(1,3,2)
plot(k_vals, maxU, '-o', 'LineWidth', 1.5);
xlabel('$k$', 'Interpreter','latex'); ylabel('\textbf{max \{u\}}', 'Interpreter','latex');
title('\textbf{Efecto de $k$ en la solucion sobre el dominio tipo $L$}', 'Interpreter','latex');
grid on; hold on; view(2); axis square;

%% Ejemplo de Dominio Tipo T

%%% --- Parámetros de la malla ---
n_x = 30; n_y = 30;
[x_v, y_v] = meshgrid(linspace(0,3,n_x+1), linspace(0,3,n_y+1));
X = x_v(:); Y = y_v(:);
nodos = [X Y];

fuera_dominio = (nodos(:,1) > 2 & nodos(:,2) < 2) | (nodos(:,1) < 1 & nodos(:,2) < 2);
nodos = nodos(~fuera_dominio,:);

tol = 1e-12;

is_gamma1 = abs(nodos(:,2)) < tol & nodos(:,1) >= 1 & nodos(:,1) <= 2;
is_gamma2 = abs(nodos(:,1)-1) < tol & nodos(:,2) >= 0 & nodos(:,2) <= 2;
is_gamma3 = abs(nodos(:,2)-2) < tol & nodos(:,1) >= 0 & nodos(:,1) <= 1;
is_gamma4 = abs(nodos(:,1)) < tol & nodos(:,2) >= 2 & nodos(:,2) <= 3;
is_gamma5 = abs(nodos(:,2)-3) < tol & nodos(:,1) >= 0 & nodos(:,1) <= 3;
is_gamma6 = abs(nodos(:,1)-3) < tol & nodos(:,2) >= 2 & nodos(:,2) <= 3;
is_gamma7 = abs(nodos(:,2)-2) < tol & nodos(:,1) >= 2 & nodos(:,1) <= 3;
is_gamma8 = abs(nodos(:,1)-2) < tol & nodos(:,2) >= 0 & nodos(:,2) <= 2;

triang = delaunay(nodos);
n_nodos = size(nodos,1);

is_edge = is_gamma2 | is_gamma3 | is_gamma7 | is_gamma8;
triang_eliminar = all(is_edge(triang), 2);
triang(triang_eliminar, :) = [];

n_triang = size(triang,1);

figure(1)
subplot(1,3,3)
trimesh(triang, nodos(:,1), nodos(:,2), zeros(size(nodos,1),1))
%shading interp; camproj perspective; axis vis3d; view(45,30)
xlabel('x'), ylabel('y'); hold on;
title('\textbf{Dominio tipo $T$}', 'Interpreter','latex')
scatter(1.5, 0.1, 'red', '*', LineWidth = 3)
view(2); axis equal tight;

%%% --- Parámetros del problema ---
k_fun = @(x,y) 1; % conductividad
c_fun = @(x,y) 0; % término de reacción
% r_fun = @(x,y) 10*exp(-50*((x-0.5).^2+(y-0.5).^2)); % fuente (como en tus códigos)
r_fun = @(x,y) 125*exp(-10*((x-1.5).^2+(y-0.1).^2)); % Fuente Gaussiana

%%% --- Selección de pared Neumann ---
% Opciones: 'gamma1','gamma2','gamma3','gamma4','gamma5','gamma6','gamma7','gamma8'
neumann_sides = ["gamma1","gamma2","gamma3","gamma5","gamma7","gamma8"];

%%% --- Ensamblaje global ---
I = []; J = []; V = []; R = zeros(n_nodos,1);

for t = 1:n_triang
id = triang(t,:);
puntos = nodos(id,:);
x1 = puntos(1,1); y1 = puntos(1,2);
x2 = puntos(2,1); y2 = puntos(2,2);
x3 = puntos(3,1); y3 = puntos(3,2);

Area = abs(det([x2-x1 x3-x1; y2-y1 y3-y1]))/2;
B = [y2-y3 y3-y1 y1-y2; x3-x2 x1-x3 x2-x1]/(2*Area);
x_m = mean(puntos(:,1)); y_m = mean(puntos(:,2));

k = k_fun(x_m, y_m);
Ke = k*Area*(B'*B);

c = c_fun(x_m, y_m);
Me = (c*Area/12)*[2 1 1; 1 2 1; 1 1 2];

rmid = r_fun(x_m, y_m);
Re = rmid*Area/3*ones(3,1);

for a = 1:3
for b = 1:3
I(end+1) = id(a);
J(end+1) = id(b);
V(end+1) = Ke(a,b) + Me(a,b);
end
R(id(a)) = R(id(a)) + Re(a);
end
end

A = sparse(I, J, V, n_nodos, n_nodos);
b = R;

%%% --- Identificación de nodos de contorno y clasificación Dirichlet/Neumann ---
tol = 1e-12;

% Por defecto: Dirichlet en todas las paredes
dirichlet_nodes = is_gamma1 | is_gamma2 | is_gamma3 |...
is_gamma4 | is_gamma5 | is_gamma6 |...
is_gamma7 | is_gamma8;
neumann_nodes = zeros(length(is_gamma1),1);

% Reemplazar la pared seleccionada por Neumann
for i = 1:numel(neumann_sides)

neumann_side = neumann_sides(i);

switch lower(neumann_side)
case 'gamma1'
dirichlet_nodes(is_gamma1) = false;
neumann_nodes = neumann_nodes + is_gamma1;
case 'gamma2'
dirichlet_nodes(is_gamma2) = false;
neumann_nodes = neumann_nodes + is_gamma2;
case 'gamma3'
dirichlet_nodes(is_gamma3) = false;
neumann_nodes = neumann_nodes + is_gamma3;
case 'gamma4'
dirichlet_nodes(is_gamma4) = false;
neumann_nodes = neumann_nodes + is_gamma4;
case 'gamma5'
dirichlet_nodes(is_gamma5) = false;
neumann_nodes = neumann_nodes + is_gamma5;
case 'gamma6'
dirichlet_nodes(is_gamma6) = false;
neumann_nodes = neumann_nodes + is_gamma6;
case 'gamma7'
dirichlet_nodes(is_gamma7) = false;
neumann_nodes = neumann_nodes + is_gamma7;
case 'gamma8'
dirichlet_nodes(is_gamma8) = false;
neumann_nodes = neumann_nodes + is_gamma8;
otherwise
error('neumann_side no reconocido. Usa gamma1,gamma2,gamma3,gamma4,gamma5,gamma6,gamma7,gamma8.')
end
end

ind_dir = find(dirichlet_nodes);
ind_neu = find(neumann_nodes);

%%% --- Aplicación de condiciones Dirichlet homogéneas (en nodos dirichlet_nodes) ---
% Igual que en tu primer código: anular filas/columnas y poner 1 en diagonal
A(ind_dir, :) = 0; A(:, ind_dir) = 0;
A(sub2ind(size(A), ind_dir, ind_dir)) = 1;
b(ind_dir) = 0;

%%% --- Resolución del sistema ---
u = A\b;
u = u(:);

%%% --- Visualización ---
figure(2)
subplot(3,3,7)
trisurf(triang, nodos(:,1), nodos(:,2), u, 'EdgeColor', 'none');
view(2); axis equal tight; colorbar;
title('\textbf{Solucion FEM (vista superior) sobre el dominio tipo $T$}', 'Interpreter','latex')
xlabel('x'), ylabel('y')

subplot(3,3,8)
trisurf(triang, nodos(:,1), nodos(:,2), u)
shading interp; camproj perspective; axis vis3d; view(125,30)
xlabel('x'), ylabel('y'), zlabel('u')
title('\textbf{Solucion FEM (superficie) sobre el dominio tipo $T$}', 'Interpreter','latex')
xlabel('x'), ylabel('y')

% Mostrar en la malla qué nodos son Dirichlet/Neumann
subplot(3,3,9)
plot(nodos(:,1), nodos(:,2), '.k'); hold on
plot(nodos(ind_dir,1), nodos(ind_dir,2), 'or', 'MarkerFaceColor','r');
plot(nodos(ind_neu,1), nodos(ind_neu,2), 'sb', 'MarkerFaceColor','b');
legend('Interior/otros','Dirichlet','Neumann','Location','bestoutside');
xlabel('x'), ylabel('y');
title('\textbf{Nodos de contorno sobre el dominio tipo $T$}', 'Interpreter','latex');
view(2); axis equal tight; xlabel('x'), ylabel('y')

%%% --- Barrido de parámetros k (igual que en tus scripts originales) ---
k_vals = [0.1 1 5 10 15 20];
maxU = zeros(size(k_vals));

for i = 1:length(k_vals)
k_fun = @(x,y) k_vals(i);

% Re-ensamblar para cada k
I = []; J = []; V = []; R = zeros(n_nodos,1);
for t = 1:n_triang
id = triang(t,:);
puntos = nodos(id,:);
x1 = puntos(1,1); y1 = puntos(1,2);
x2 = puntos(2,1); y2 = puntos(2,2);
x3 = puntos(3,1); y3 = puntos(3,2);

Area = abs(det([x2-x1 x3-x1; y2-y1 y3-y1]))/2;
B = [y2-y3 y3-y1 y1-y2; x3-x2 x1-x3 x2-x1]/(2*Area);
x_m = mean(puntos(:,1)); y_m = mean(puntos(:,2));

k = k_fun(x_m, y_m);
Ke = k*Area*(B'*B);
c = c_fun(x_m, y_m);
Me = (c*Area/12)*[2 1 1; 1 2 1; 1 1 2];
rmid = r_fun(x_m, y_m);
Re = rmid*Area/3*ones(3,1);

for a = 1:3
for b = 1:3
I(end+1) = id(a);
J(end+1) = id(b);
V(end+1) = Ke(a,b) + Me(a,b);
end
R(id(a)) = R(id(a)) + Re(a);
end
end
A_k = sparse(I, J, V, n_nodos, n_nodos);
b_k = R;

% Impón Dirichlet (mismo conjunto de nodos)
A_k(ind_dir, :) = 0; A_k(:, ind_dir) = 0;
A_k(sub2ind(size(A_k), ind_dir, ind_dir)) = 1;
b_k(ind_dir) = 0;

% Resolver y almacenar max
u_k = A_k \ b_k;
maxU(i) = max(u_k);
end

figure(3)
subplot(1,3,3)
plot(k_vals, maxU, '-o', 'LineWidth', 1.5);
xlabel('$k$', 'Interpreter','latex'); ylabel('\textbf{max \{u\}}', 'Interpreter','latex');
title('\textbf{Efecto de $k$ en la solucion sobre el dominio tipo $T$}', 'Interpreter','latex');
grid on; hold on; view(2); axis square;
</source>

[[Categoría:MNEDP|MNEDP]]
[[Categoría:MNEDP25/26|2025-26]]

Archivo:Poster MNEDP DIF.jpg

2025-11-12T10:09:06Z

Daniel Rodríguez Calderón: Trabajo de Métodos Numéricos en Ecuaciones Diferenciales en Derivadas Parciales sobre la Ecuación del Calor Estacionaria en formato póster realizado por el grupo DIF.

Trabajo de Métodos Numéricos en Ecuaciones Diferenciales en Derivadas Parciales sobre la Ecuación del Calor Estacionaria en formato póster realizado por el grupo DIF.

Ecuación del Calor Estacionaria (Grupo DIF)

2025-11-11T09:04:30Z

Daniel Rodríguez Calderón:

{{ TrabajoED | Ecuación del Calor Estacionaria (Grupo DIF) | [[:Categoría:MNEDP|MNEDP]]|[[:Categoría:MNEDP25/26|2025-26]] |
*Daniel Rodríguez Calderón,
*Israel López Moreno,
*Francisco Lavao Amengual.}}

=Póster del Trabajo=

{|style="margin: 0 auto;"
| [[Archivo:Poster_Ecuacion_Poisson_ADMR.jpeg|900px|thumb|right|Trabajo de Métodos Numéricos en Ecuaciones Diferenciales en Derivadas Parciales sobre la Ecuación del Calor Estacionaria en formato póster realizado por el grupo DIF.]]
|}

=Códigos de Matlab=
<source lang="matlab">
%% MEF de la Ecuación del Calor con condiciones Dirichlet homogéneas + Neumann homogéneas
clc; clear all; close all

%% Ejemplo de Dominio Cuadrado Unidad

%%% --- Parámetros de la malla ---
n_x = 30; n_y = 30; tol = 1e-12;

[x_v, y_v] = meshgrid(linspace(0,1,n_x+1), linspace(0,1,n_y+1));
X = x_v(:); Y = y_v(:);
nodos = [X Y];

triang = delaunay(nodos);
n_nodos = size(nodos,1);
n_triang = size(triang,1);

figure(1)
subplot(1,3,1)
trimesh(triang, nodos(:,1), nodos(:,2), zeros(size(nodos,1),1))
hold on; xlabel('x'), ylabel('y');
%shading interp; camproj perspective; axis vis3d; view(45,30)
title('\textbf{Dominio cuadrado unidad}', 'Interpreter','latex')
scatter(0.1, 0.5, 'red', '*', LineWidth = 3)
view(2); axis equal tight;

%%% --- Parámetros del problema ---
k_fun = @(x,y) 1; % conductividad
c_fun = @(x,y) 0; % término de reacción
% r_fun = @(x,y) 10*exp(-50*((x-0.5).^2+(y-0.5).^2)); % fuente (como en tus códigos)
r_fun = @(x,y) 350*exp(-10*((x-0.1).^2+(y-0.5).^2)); % Fuente Gaussiana

%%% --- Selección de frontera Neumann ---
% Opciones: 'gamma1','gamma2','gamma3','gamma4'
neumann_sides = ["gamma1","gamma2","gamma3"];
g_val = 0; % flujo de Neumann en la pared seleccionada (0 según tu petición)

%%% --- Ensamblaje global ---
I = []; J = []; V = []; R = zeros(n_nodos,1);

for t = 1:n_triang
id = triang(t,:);
puntos = nodos(id,:);
x1 = puntos(1,1); y1 = puntos(1,2);
x2 = puntos(2,1); y2 = puntos(2,2);
x3 = puntos(3,1); y3 = puntos(3,2);

Area = abs(det([x2-x1 x3-x1; y2-y1 y3-y1]))/2;
B = [y2-y3 y3-y1 y1-y2; x3-x2 x1-x3 x2-x1]/(2*Area);
x_m = mean(puntos(:,1)); y_m = mean(puntos(:,2));

k = k_fun(x_m, y_m);
Ke = k*Area*(B'*B);

c = c_fun(x_m, y_m);
Me = (c*Area/12)*[2 1 1; 1 2 1; 1 1 2];

rmid = r_fun(x_m, y_m);
Re = rmid*Area/3*ones(3,1);

for a = 1:3
for b = 1:3
I(end+1) = id(a);
J(end+1) = id(b);
V(end+1) = Ke(a,b) + Me(a,b);
end
R(id(a)) = R(id(a)) + Re(a);
end
end

A = sparse(I, J, V, n_nodos, n_nodos);
b = R;

%%% --- Identificación de nodos de contorno y clasificación Dirichlet/Neumann ---
is_gamma1 = abs(nodos(:,2)) < tol & nodos(:,1) >= 0 & nodos(:,1) <= 1;
is_gamma2 = abs(nodos(:,1)) < tol & nodos(:,2) >= 0 & nodos(:,2) <= 1;
is_gamma3 = abs(nodos(:,2)-1) < tol & nodos(:,1) >= 0 & nodos(:,1) <= 1;
is_gamma4 = abs(nodos(:,1)-1) < tol & nodos(:,2) >= 0 & nodos(:,2) <= 1;

% Por defecto: Dirichlet en todas las paredes
dirichlet_nodes = is_gamma1 | is_gamma2 | is_gamma3 | is_gamma4;
neumann_nodes = zeros(length(is_gamma1), 1);

% Reemplazar la pared seleccionada por Neumann
for i = 1:numel(neumann_sides)
neumann_side = neumann_sides(i);

switch lower(neumann_side)
case 'gamma1'
dirichlet_nodes(is_gamma1) = false;
neumann_nodes = neumann_nodes + is_gamma1;
case 'gamma2'
dirichlet_nodes(is_gamma2) = false;
neumann_nodes = neumann_nodes + is_gamma2;
case 'gamma3'
dirichlet_nodes(is_gamma3) = false;
neumann_nodes = neumann_nodes + is_gamma3;
case 'gamma4'
dirichlet_nodes(is_gamma4) = false;
neumann_nodes = neumann_nodes + is_gamma4;
otherwise
error('neumann_side no reconocido. Usa gamma1,gamma2,gamma3,gamma4,gamma5,gamma6.')
end
end

ind_dir = find(dirichlet_nodes);
ind_neu = find(neumann_nodes);

%%% --- Aplicación de condiciones Dirichlet homogéneas (en nodos dirichlet_nodes) ---
% Igual que en tu primer código: anular filas/columnas y poner 1 en diagonal
A(ind_dir, :) = 0; A(:, ind_dir) = 0;
A(sub2ind(size(A), ind_dir, ind_dir)) = 1;
b(ind_dir) = 0;

%%% --- Resolución del sistema ---
u = A\b;
u = u(:);

%%% --- Visualización ---
figure(2);
subplot(3,3,1)
trisurf(triang, nodos(:,1), nodos(:,2), u, 'EdgeColor', 'none');
view(2); axis equal tight; colorbar; colormap(turbo);
title('\textbf{Solucion FEM (vista superior) sobre el dominio cuadrado unidad}', 'Interpreter','latex')
xlabel('x'), ylabel('y')

subplot(3,3,2)
trisurf(triang, nodos(:,1), nodos(:,2), u)
shading interp; camproj perspective; axis vis3d; view(45,30)
xlabel('x'), ylabel('y'), zlabel('u')
title('\textbf{Solucion FEM (superficie) sobre el dominio cuadrado unidad}', 'Interpreter','latex')

% Mostrar en la malla qué nodos son Dirichlet/Neumann
subplot(3,3,3)
plot(nodos(:,1), nodos(:,2), '.k'); hold on
plot(nodos(ind_dir,1), nodos(ind_dir,2), 'or', 'MarkerFaceColor','r');
plot(nodos(ind_neu,1), nodos(ind_neu,2), 'sb', 'MarkerFaceColor','b');
legend('Interior/otros','Dirichlet','Neumann','Location','bestoutside')
view(2); axis equal tight; xlabel('x'), ylabel('y');
title('\textbf{Nodos de contorno sobre el dominio cuadrado unidad}', 'Interpreter','latex');

%%% --- Barrido de parámetros k (igual que en tus scripts originales) ---
k_vals = [0.1 1 5 10 15 20];
maxU = zeros(size(k_vals));

for i = 1:length(k_vals)
k_fun = @(x,y) k_vals(i);

% Re-ensamblar para cada k
I = []; J = []; V = []; R = zeros(n_nodos,1);
for t = 1:n_triang
id = triang(t,:);
puntos = nodos(id,:);
x1 = puntos(1,1); y1 = puntos(1,2);
x2 = puntos(2,1); y2 = puntos(2,2);
x3 = puntos(3,1); y3 = puntos(3,2);

Area = abs(det([x2-x1 x3-x1; y2-y1 y3-y1]))/2;
B = [y2-y3 y3-y1 y1-y2; x3-x2 x1-x3 x2-x1]/(2*Area);
x_m = mean(puntos(:,1)); y_m = mean(puntos(:,2));

k = k_fun(x_m, y_m);
Ke = k*Area*(B'*B);
c = c_fun(x_m, y_m);
Me = (c*Area/12)*[2 1 1; 1 2 1; 1 1 2];
rmid = r_fun(x_m, y_m);
Re = rmid*Area/3*ones(3,1);

for a = 1:3
for b = 1:3
I(end+1) = id(a);
J(end+1) = id(b);
V(end+1) = Ke(a,b) + Me(a,b);
end
R(id(a)) = R(id(a)) + Re(a);
end
end
A_k = sparse(I, J, V, n_nodos, n_nodos);
b_k = R;

% Impón Dirichlet (mismo conjunto de nodos)
A_k(ind_dir, :) = 0; A_k(:, ind_dir) = 0;
A_k(sub2ind(size(A_k), ind_dir, ind_dir)) = 1;
b_k(ind_dir) = 0;

% Resolver y almacenar max
u_k = A_k \ b_k;
maxU(i) = max(u_k);
end

figure(3)
subplot(1,3,1)
plot(k_vals, maxU, '-o', 'LineWidth', 1.5);
xlabel('$k$', 'Interpreter','latex'); ylabel('\textbf{max \{u\}}', 'Interpreter','latex');
title('\textbf{Efecto de $k$ en la solucion sobre el dominio cuadrado unidad}', 'Interpreter','latex');
grid on; hold on; view(2); axis square;

%% Ejemplo de Dominio Tipo L

%%% --- Parámetros de la malla ---
n_x = 30; n_y = 30;
[x_v, y_v] = meshgrid(linspace(0,3,n_x+1), linspace(0,2,n_y+1));
X = x_v(:); Y = y_v(:);
nodos = [X Y];

fuera_dominio = nodos(:,1) > 1 & nodos(:,2) < 1;
nodos = nodos(~fuera_dominio,:);

tol = 1e-12;

is_gamma1 = abs(nodos(:,2)) < tol & nodos(:,1) >= 0 & nodos(:,1) <= 1;
is_gamma2 = abs(nodos(:,1)) < tol & nodos(:,2) >= 0 & nodos(:,2) <= 2;
is_gamma3 = abs(nodos(:,2)-2) < tol & nodos(:,1) >= 0 & nodos(:,1) <= 3;
is_gamma4 = abs(nodos(:,1)- 3) < tol & nodos(:,2) >= 1 & nodos(:,2) <= 2;
is_gamma5 = abs(nodos(:,2)-1) < tol & nodos(:,1) >= 1 & nodos(:,1) <= 3;
is_gamma6 = abs(nodos(:,1)-1) < tol & nodos(:,2) >= 0 & nodos(:,2) <= 1;
is_beta = nodos(:,1) > 1 & nodos(:,1) <= 3 & nodos(:,2) >= 0 & nodos(:,2) < 1;

triang = delaunay(nodos);
n_nodos = size(nodos,1);

is_edge = is_gamma5 | is_gamma6;
triang_eliminar = all(is_edge(triang), 2);
triang(triang_eliminar, :) = [];

n_triang = size(triang,1);

figure(1)
subplot(1,3,2)
trimesh(triang, nodos(:,1), nodos(:,2), zeros(size(nodos,1),1))
%shading interp; camproj perspective; axis vis3d; view(45,30)
xlabel('x'), ylabel('y'); hold on;
title('\textbf{Dominio tipo $L$}', 'Interpreter','latex')
scatter([0.5 1.5], [0.1 1.9], 'red', '*', LineWidth = 3)
view(2); axis equal tight;

%%% --- Parámetros del problema ---
k_fun = @(x,y) 1; % conductividad
c_fun = @(x,y) 0; % término de reacción
% r_fun = @(x,y) 10*exp(-50*((x-0.5).^2+(y-0.5).^2)); % fuente (como en tus códigos)
r_fun = @(x,y) 75*exp(-10*((x-0.5).^2+(y-0.1).^2)) + 75*exp(-10*((x-1.5).^2+(y-1.9).^2)); % Fuente Gaussiana

%%% --- Selección de frontera Neumann ---
% Opciones: 'gamma1','gamma2','gamma3','gamma4','gamma5','gamma6'
neumann_sides = ["gamma1","gamma2","gamma3","gamma5","gamma6"];
g_val = 0; % flujo de Neumann en la pared seleccionada (0 según tu petición)

%%% --- Ensamblaje global ---
I = []; J = []; V = []; R = zeros(n_nodos,1);

for t = 1:n_triang
id = triang(t,:);
puntos = nodos(id,:);
x1 = puntos(1,1); y1 = puntos(1,2);
x2 = puntos(2,1); y2 = puntos(2,2);
x3 = puntos(3,1); y3 = puntos(3,2);

Area = abs(det([x2-x1 x3-x1; y2-y1 y3-y1]))/2;
B = [y2-y3 y3-y1 y1-y2; x3-x2 x1-x3 x2-x1]/(2*Area);
x_m = mean(puntos(:,1)); y_m = mean(puntos(:,2));

k = k_fun(x_m, y_m);
Ke = k*Area*(B'*B);

c = c_fun(x_m, y_m);
Me = (c*Area/12)*[2 1 1; 1 2 1; 1 1 2];

rmid = r_fun(x_m, y_m);
Re = rmid*Area/3*ones(3,1);

for a = 1:3
for b = 1:3
I(end+1) = id(a);
J(end+1) = id(b);
V(end+1) = Ke(a,b) + Me(a,b);
end
R(id(a)) = R(id(a)) + Re(a);
end
end

A = sparse(I, J, V, n_nodos, n_nodos);
b = R;

%%% --- Identificación de nodos de contorno y clasificación Dirichlet/Neumann ---
is_gamma1 = abs(nodos(:,2)) < tol & nodos(:,1) >= 0 & nodos(:,1) <= 1;
is_gamma2 = abs(nodos(:,1)) < tol & nodos(:,2) >= 0 & nodos(:,2) <= 2;
is_gamma3 = abs(nodos(:,2)-2) < tol & nodos(:,1) >= 0 & nodos(:,1) <= 3;
is_gamma4 = abs(nodos(:,1)- 3) < tol & nodos(:,2) >= 1 & nodos(:,2) <= 2;
is_gamma5 = abs(nodos(:,2)-1) < tol & nodos(:,1) >= 1 & nodos(:,1) <= 3;
is_gamma6 = abs(nodos(:,1)-1) < tol & nodos(:,2) >= 0 & nodos(:,2) <= 1;
is_beta = nodos(:,1) > 1 & nodos(:,1) <= 3 & nodos(:,2) >= 0 & nodos(:,2) < 1;

% Por defecto: Dirichlet en todas las paredes
dirichlet_nodes = is_gamma1 | is_gamma2 | is_gamma3 |...
is_gamma4 | is_gamma5 | is_gamma6 | is_beta;
neumann_nodes = zeros(length(is_gamma1),1);

% Reemplazar la pared seleccionada por Neumann
for i = 1:numel(neumann_sides)
neumann_side = neumann_sides(i);

switch lower(neumann_side)
case 'gamma1'
dirichlet_nodes(is_gamma1) = false;
neumann_nodes = neumann_nodes + is_gamma1;
case 'gamma2'
dirichlet_nodes(is_gamma2) = false;
neumann_nodes = neumann_nodes + is_gamma2;
case 'gamma3'
dirichlet_nodes(is_gamma3) = false;
neumann_nodes = neumann_nodes + is_gamma3;
case 'gamma4'
dirichlet_nodes(is_gamma4) = false;
neumann_nodes = neumann_nodes + is_gamma4;
case 'gamma5'
dirichlet_nodes(is_gamma5) = false;
neumann_nodes = neumann_nodes + is_gamma5;
case 'gamma6'
dirichlet_nodes(is_gamma6) = false;
neumann_nodes = neumann_nodes + is_gamma6;
otherwise
error('neumann_side no reconocido. Usa gamma1,gamma2,gamma3,gamma4,gamma5,gamma6.')
end
end

ind_dir = find(dirichlet_nodes);
ind_neu = find(neumann_nodes);

%%% --- Aplicación de condiciones Dirichlet homogéneas (en nodos dirichlet_nodes) ---
% Igual que en tu primer código: anular filas/columnas y poner 1 en diagonal
A(ind_dir, :) = 0; A(:, ind_dir) = 0;
A(sub2ind(size(A), ind_dir, ind_dir)) = 1;
b(ind_dir) = 0;

%%% --- Resolución del sistema ---
u = A\b;
u = u(:);

%%% --- Visualización ---
figure(2);
subplot(3,3,4)
trisurf(triang, nodos(:,1), nodos(:,2), u, 'EdgeColor', 'none');
view(2); axis equal tight; colorbar;
title('\textbf{Solucion FEM (vista superior) sobre el dominio tipo $L$}', 'Interpreter','latex')
xlabel('x'), ylabel('y')

subplot(3,3,5)
trisurf(triang, nodos(:,1), nodos(:,2), u)
shading interp; camproj perspective; axis vis3d; view(45,30)
xlabel('x'), ylabel('y'), zlabel('u')
title('\textbf{Solucion FEM (superficie) sobre el dominio tipo $L$}', 'Interpreter','latex')

% Mostrar en la malla qué nodos son Dirichlet/Neumann
subplot(3,3,6)
plot(nodos(:,1), nodos(:,2), '.k'); hold on
plot(nodos(ind_dir,1), nodos(ind_dir,2), 'or', 'MarkerFaceColor','r');
plot(nodos(ind_neu,1), nodos(ind_neu,2), 'sb', 'MarkerFaceColor','b');
legend('Interior/otros','Dirichlet','Neumann','Location','bestoutside')
view(2); axis equal tight; xlabel('x'), ylabel('y')
title('\textbf{Nodos de contorno sobre el dominio tipo $L$}', 'Interpreter','latex');

%%% --- Barrido de parámetros k (igual que en tus scripts originales) ---
k_vals = [0.1 1 5 10 15 20];
maxU = zeros(size(k_vals));

for i = 1:length(k_vals)
k_fun = @(x,y) k_vals(i);

% Re-ensamblar para cada k
I = []; J = []; V = []; R = zeros(n_nodos,1);
for t = 1:n_triang
id = triang(t,:);
puntos = nodos(id,:);
x1 = puntos(1,1); y1 = puntos(1,2);
x2 = puntos(2,1); y2 = puntos(2,2);
x3 = puntos(3,1); y3 = puntos(3,2);

Area = abs(det([x2-x1 x3-x1; y2-y1 y3-y1]))/2;
B = [y2-y3 y3-y1 y1-y2; x3-x2 x1-x3 x2-x1]/(2*Area);
x_m = mean(puntos(:,1)); y_m = mean(puntos(:,2));

k = k_fun(x_m, y_m);
Ke = k*Area*(B'*B);
c = c_fun(x_m, y_m);
Me = (c*Area/12)*[2 1 1; 1 2 1; 1 1 2];
rmid = r_fun(x_m, y_m);
Re = rmid*Area/3*ones(3,1);

for a = 1:3
for b = 1:3
I(end+1) = id(a);
J(end+1) = id(b);
V(end+1) = Ke(a,b) + Me(a,b);
end
R(id(a)) = R(id(a)) + Re(a);
end
end
A_k = sparse(I, J, V, n_nodos, n_nodos);
b_k = R;

% Impón Dirichlet (mismo conjunto de nodos)
A_k(ind_dir, :) = 0; A_k(:, ind_dir) = 0;
A_k(sub2ind(size(A_k), ind_dir, ind_dir)) = 1;
b_k(ind_dir) = 0;

% Resolver y almacenar max
u_k = A_k \ b_k;
maxU(i) = max(u_k);
end

figure(3)
subplot(1,3,2)
plot(k_vals, maxU, '-o', 'LineWidth', 1.5);
xlabel('$k$', 'Interpreter','latex'); ylabel('\textbf{max \{u\}}', 'Interpreter','latex');
title('\textbf{Efecto de $k$ en la solucion sobre el dominio tipo $L$}', 'Interpreter','latex');
grid on; hold on; view(2); axis square;

%% Ejemplo de Dominio Tipo T

%%% --- Parámetros de la malla ---
n_x = 30; n_y = 30;
[x_v, y_v] = meshgrid(linspace(0,3,n_x+1), linspace(0,3,n_y+1));
X = x_v(:); Y = y_v(:);
nodos = [X Y];

fuera_dominio = (nodos(:,1) > 2 & nodos(:,2) < 2) | (nodos(:,1) < 1 & nodos(:,2) < 2);
nodos = nodos(~fuera_dominio,:);

tol = 1e-12;

is_gamma1 = abs(nodos(:,2)) < tol & nodos(:,1) >= 1 & nodos(:,1) <= 2;
is_gamma2 = abs(nodos(:,1)-1) < tol & nodos(:,2) >= 0 & nodos(:,2) <= 2;
is_gamma3 = abs(nodos(:,2)-2) < tol & nodos(:,1) >= 0 & nodos(:,1) <= 1;
is_gamma4 = abs(nodos(:,1)) < tol & nodos(:,2) >= 2 & nodos(:,2) <= 3;
is_gamma5 = abs(nodos(:,2)-3) < tol & nodos(:,1) >= 0 & nodos(:,1) <= 3;
is_gamma6 = abs(nodos(:,1)-3) < tol & nodos(:,2) >= 2 & nodos(:,2) <= 3;
is_gamma7 = abs(nodos(:,2)-2) < tol & nodos(:,1) >= 2 & nodos(:,1) <= 3;
is_gamma8 = abs(nodos(:,1)-2) < tol & nodos(:,2) >= 0 & nodos(:,2) <= 2;

triang = delaunay(nodos);
n_nodos = size(nodos,1);

is_edge = is_gamma2 | is_gamma3 | is_gamma7 | is_gamma8;
triang_eliminar = all(is_edge(triang), 2);
triang(triang_eliminar, :) = [];

n_triang = size(triang,1);

figure(1)
subplot(1,3,3)
trimesh(triang, nodos(:,1), nodos(:,2), zeros(size(nodos,1),1))
%shading interp; camproj perspective; axis vis3d; view(45,30)
xlabel('x'), ylabel('y'); hold on;
title('\textbf{Dominio tipo $T$}', 'Interpreter','latex')
scatter(1.5, 0.1, 'red', '*', LineWidth = 3)
view(2); axis equal tight;

%%% --- Parámetros del problema ---
k_fun = @(x,y) 1; % conductividad
c_fun = @(x,y) 0; % término de reacción
% r_fun = @(x,y) 10*exp(-50*((x-0.5).^2+(y-0.5).^2)); % fuente (como en tus códigos)
r_fun = @(x,y) 125*exp(-10*((x-1.5).^2+(y-0.1).^2)); % Fuente Gaussiana

%%% --- Selección de pared Neumann ---
% Opciones: 'gamma1','gamma2','gamma3','gamma4','gamma5','gamma6','gamma7','gamma8'
neumann_sides = ["gamma1","gamma2","gamma3","gamma5","gamma7","gamma8"];

%%% --- Ensamblaje global ---
I = []; J = []; V = []; R = zeros(n_nodos,1);

for t = 1:n_triang
id = triang(t,:);
puntos = nodos(id,:);
x1 = puntos(1,1); y1 = puntos(1,2);
x2 = puntos(2,1); y2 = puntos(2,2);
x3 = puntos(3,1); y3 = puntos(3,2);

Area = abs(det([x2-x1 x3-x1; y2-y1 y3-y1]))/2;
B = [y2-y3 y3-y1 y1-y2; x3-x2 x1-x3 x2-x1]/(2*Area);
x_m = mean(puntos(:,1)); y_m = mean(puntos(:,2));

k = k_fun(x_m, y_m);
Ke = k*Area*(B'*B);

c = c_fun(x_m, y_m);
Me = (c*Area/12)*[2 1 1; 1 2 1; 1 1 2];

rmid = r_fun(x_m, y_m);
Re = rmid*Area/3*ones(3,1);

for a = 1:3
for b = 1:3
I(end+1) = id(a);
J(end+1) = id(b);
V(end+1) = Ke(a,b) + Me(a,b);
end
R(id(a)) = R(id(a)) + Re(a);
end
end

A = sparse(I, J, V, n_nodos, n_nodos);
b = R;

%%% --- Identificación de nodos de contorno y clasificación Dirichlet/Neumann ---
tol = 1e-12;

% Por defecto: Dirichlet en todas las paredes
dirichlet_nodes = is_gamma1 | is_gamma2 | is_gamma3 |...
is_gamma4 | is_gamma5 | is_gamma6 |...
is_gamma7 | is_gamma8;
neumann_nodes = zeros(length(is_gamma1),1);

% Reemplazar la pared seleccionada por Neumann
for i = 1:numel(neumann_sides)

neumann_side = neumann_sides(i);

switch lower(neumann_side)
case 'gamma1'
dirichlet_nodes(is_gamma1) = false;
neumann_nodes = neumann_nodes + is_gamma1;
case 'gamma2'
dirichlet_nodes(is_gamma2) = false;
neumann_nodes = neumann_nodes + is_gamma2;
case 'gamma3'
dirichlet_nodes(is_gamma3) = false;
neumann_nodes = neumann_nodes + is_gamma3;
case 'gamma4'
dirichlet_nodes(is_gamma4) = false;
neumann_nodes = neumann_nodes + is_gamma4;
case 'gamma5'
dirichlet_nodes(is_gamma5) = false;
neumann_nodes = neumann_nodes + is_gamma5;
case 'gamma6'
dirichlet_nodes(is_gamma6) = false;
neumann_nodes = neumann_nodes + is_gamma6;
case 'gamma7'
dirichlet_nodes(is_gamma7) = false;
neumann_nodes = neumann_nodes + is_gamma7;
case 'gamma8'
dirichlet_nodes(is_gamma8) = false;
neumann_nodes = neumann_nodes + is_gamma8;
otherwise
error('neumann_side no reconocido. Usa gamma1,gamma2,gamma3,gamma4,gamma5,gamma6,gamma7,gamma8.')
end
end

ind_dir = find(dirichlet_nodes);
ind_neu = find(neumann_nodes);

%%% --- Aplicación de condiciones Dirichlet homogéneas (en nodos dirichlet_nodes) ---
% Igual que en tu primer código: anular filas/columnas y poner 1 en diagonal
A(ind_dir, :) = 0; A(:, ind_dir) = 0;
A(sub2ind(size(A), ind_dir, ind_dir)) = 1;
b(ind_dir) = 0;

%%% --- Resolución del sistema ---
u = A\b;
u = u(:);

%%% --- Visualización ---
figure(2)
subplot(3,3,7)
trisurf(triang, nodos(:,1), nodos(:,2), u, 'EdgeColor', 'none');
view(2); axis equal tight; colorbar;
title('\textbf{Solucion FEM (vista superior) sobre el dominio tipo $T$}', 'Interpreter','latex')
xlabel('x'), ylabel('y')

subplot(3,3,8)
trisurf(triang, nodos(:,1), nodos(:,2), u)
shading interp; camproj perspective; axis vis3d; view(125,30)
xlabel('x'), ylabel('y'), zlabel('u')
title('\textbf{Solucion FEM (superficie) sobre el dominio tipo $T$}', 'Interpreter','latex')
xlabel('x'), ylabel('y')

% Mostrar en la malla qué nodos son Dirichlet/Neumann
subplot(3,3,9)
plot(nodos(:,1), nodos(:,2), '.k'); hold on
plot(nodos(ind_dir,1), nodos(ind_dir,2), 'or', 'MarkerFaceColor','r');
plot(nodos(ind_neu,1), nodos(ind_neu,2), 'sb', 'MarkerFaceColor','b');
legend('Interior/otros','Dirichlet','Neumann','Location','bestoutside');
xlabel('x'), ylabel('y');
title('\textbf{Nodos de contorno sobre el dominio tipo $T$}', 'Interpreter','latex');
view(2); axis equal tight; xlabel('x'), ylabel('y')

%%% --- Barrido de parámetros k (igual que en tus scripts originales) ---
k_vals = [0.1 1 5 10 15 20];
maxU = zeros(size(k_vals));

for i = 1:length(k_vals)
k_fun = @(x,y) k_vals(i);

% Re-ensamblar para cada k
I = []; J = []; V = []; R = zeros(n_nodos,1);
for t = 1:n_triang
id = triang(t,:);
puntos = nodos(id,:);
x1 = puntos(1,1); y1 = puntos(1,2);
x2 = puntos(2,1); y2 = puntos(2,2);
x3 = puntos(3,1); y3 = puntos(3,2);

Area = abs(det([x2-x1 x3-x1; y2-y1 y3-y1]))/2;
B = [y2-y3 y3-y1 y1-y2; x3-x2 x1-x3 x2-x1]/(2*Area);
x_m = mean(puntos(:,1)); y_m = mean(puntos(:,2));

k = k_fun(x_m, y_m);
Ke = k*Area*(B'*B);
c = c_fun(x_m, y_m);
Me = (c*Area/12)*[2 1 1; 1 2 1; 1 1 2];
rmid = r_fun(x_m, y_m);
Re = rmid*Area/3*ones(3,1);

for a = 1:3
for b = 1:3
I(end+1) = id(a);
J(end+1) = id(b);
V(end+1) = Ke(a,b) + Me(a,b);
end
R(id(a)) = R(id(a)) + Re(a);
end
end
A_k = sparse(I, J, V, n_nodos, n_nodos);
b_k = R;

% Impón Dirichlet (mismo conjunto de nodos)
A_k(ind_dir, :) = 0; A_k(:, ind_dir) = 0;
A_k(sub2ind(size(A_k), ind_dir, ind_dir)) = 1;
b_k(ind_dir) = 0;

% Resolver y almacenar max
u_k = A_k \ b_k;
maxU(i) = max(u_k);
end

figure(3)
subplot(1,3,3)
plot(k_vals, maxU, '-o', 'LineWidth', 1.5);
xlabel('$k$', 'Interpreter','latex'); ylabel('\textbf{max \{u\}}', 'Interpreter','latex');
title('\textbf{Efecto de $k$ en la solucion sobre el dominio tipo $T$}', 'Interpreter','latex');
grid on; hold on; view(2); axis square;
</source>

[[Categoría:MNEDP|MNEDP]]
[[Categoría:MNEDP25/26|2025-26]]

Ecuación del Calor Estacionaria (Grupo DIF)

2025-11-10T16:14:17Z

Daniel Rodríguez Calderón:

Ecuación del Calor Estacionaria (Grupo DIF)

2025-11-10T16:13:54Z

Daniel Rodríguez Calderón:

Ecuación del Calor Estacionaria (Grupo DIF)

2025-11-10T16:08:47Z

Daniel Rodríguez Calderón:

Ecuación del Calor Estacionaria (Grupo DIF)

2025-11-10T16:03:07Z

Daniel Rodríguez Calderón:

Ecuación del Calor Estacionaria (Grupo DIF)

2025-11-10T16:02:37Z

Daniel Rodríguez Calderón:

{{ TrabajoED | Series de Fourier. Grupo 6-A | [[:Categoría:MNEDP|MNEDP]]|[[:Categoría:MNEDP25/26|2025-26]] |
*Daniel Rodríguez Calderón,
*Israel López Moreno,
*Francisco Lavao Amengual}}

[[Categoría:MNEDP|MNEDP]]
[[Categoría:MNEDP25/26|2025-26]]

Ecuación del Calor Estacionaria (Grupo DIF)

2025-11-10T16:00:43Z

Daniel Rodríguez Calderón:

{{ TrabajoED | Series de Fourier. Grupo 6-A | [[:Categoría:MNEDP|MNEDP]]|[[:Categoría:MNEDP25/26|2025-26]] |
*Daniel Rodríguez Calderón,
*Israel López Moreno,
*Francisco Lavao Amengual}}

Ecuación del Calor Estacionaria (Grupo DIF)

2025-11-10T15:59:43Z

Daniel Rodríguez Calderón: Página creada con «{{ TrabajoED | Ecuación del Calor Estacionaria (Grupo DIF) | MNEDP|2025-26 | *Daniel Rodríguez Calderón, *Israel Lópe...»

Ecuación de Laplace y de Poisson (ADMR)

2025-04-20T11:50:33Z

Daniel Rodríguez Calderón:

Ecuación de Laplace y de Poisson (ADMR)

2025-04-20T11:50:16Z

Daniel Rodríguez Calderón:

Archivo:Poster Ecuacion Poisson ADMR.jpeg

2025-04-20T11:47:36Z

Daniel Rodríguez Calderón:

Ecuación de Laplace y de Poisson (ADMR)

2025-04-20T11:40:44Z

Daniel Rodríguez Calderón:

{{ TrabajoED | Ecuación de Laplace y de Poisson (Grupo ADMR). | [[:Categoría:EDP|EDP]]|[[:Categoría:EDP24/25|2024-25]] |
*Ángel De Lucas Miranda,
*Daniel Rodríguez Calderón,
*Marcos Cabellos Hernández,
*Rafael Pascual Ortega.}}

[[Archivo:GrupoADMR Calor Errores1.jpg|750px|thumb|left|Representación gráfica de las soluciones del problema acotado (<math>a = b = 20</math>) y no acotado truncando las series de Fourier en su término <math>100</math>, para un tiempo máximo de <math>1</math>; y de sus errores. Nótese el orden de magnitud de los errores -pequeño-, los picos generados por los errores numéricos dados por la fórmula del trapecio y la suavidad en la frontera dada por las condiciones en esta.]]

[[Categoría:EDP]]
[[Categoría:EDP24/25]]

Ecuación de Laplace y de Poisson (ADMR)

2025-04-20T11:39:47Z

Daniel Rodríguez Calderón:

{{ TrabajoED | Ecuación de Laplace y de Poisson (Grupo ADMR). | [[:Categoría:EDP|EDP]]|[[:Categoría:EDP24/25|2024-25]] |
*Ángel De Lucas Miranda,
*Daniel Rodríguez Calderón,
*Marcos Cabellos Hernández,
*Rafael Pascual Ortega.}}

[[Archivo:GrupoADMR Calor Errores1.jpg|750px|thumb|left Representación gráfica de las soluciones del problema acotado (<math>a = b = 20</math>) y no acotado truncando las series de Fourier en su término <math>100</math>, para un tiempo máximo de <math>1</math>; y de sus errores. Nótese el orden de magnitud de los errores -pequeño-, los picos generados por los errores numéricos dados por la fórmula del trapecio y la suavidad en la frontera dada por las condiciones en esta.]]

{|style="margin: 0 auto;"
| [[Archivo:GrupoADMR Calor Errores DNe.jpeg|550px|thumb|upright|Podemos ver además que ambos sistemas se asemejan cada vez más, por lo que ambas aproximaciones son igualmente útiles para el caso no acotado.]]
|}

[[Categoría:EDP]]
[[Categoría:EDP24/25]]

Ecuación de Laplace y de Poisson (ADMR)

2025-04-20T11:36:59Z

Daniel Rodríguez Calderón: Página creada con «{{ TrabajoED | Ecuación de Laplace y de Poisson (Grupo ADMR). | EDP|2024-25 | *Ángel De Lucas Miranda, *Daniel Rodríguez C...»

Ecuación del calor (ADMR)

2025-03-19T16:30:27Z

Daniel Rodríguez Calderón: /* Coeficiente de difusión */

{{ TrabajoED | Ecuación del calor (Grupo ADMR). | [[:Categoría:EDP|EDP]]|[[:Categoría:EDP24/25|2024-25]] |
*Ángel De Lucas Miranda,
*Daniel Rodríguez Calderón,
*Marcos Cabellos Hernández,
*Rafael Pascual Ortega.}}

=Introducción=

El desarrollo de la ecuación del calor se atribuye principalmente al matemático y físico francés Joseph Fourier, quien en 1822 publicó su obra ''Théorie analytique de la chaleur''. En este trabajo formuló una ecuación diferencial en derivadas parciales que describe la propagación del calor por difusión en un medio continuo. En esta ocasión, profundicemos en dos aspectos de este modelo.

Desde un punto de vista físico, estas ecuaciones permiten la modelización de fenómenos como el calentamiento de materiales o la difusión de solutos en fluidos estacionarios. Sin embargo, la influencia de la física es más profunda: de hecho, el coeficiente de difusión lleva su nombre por la interpretación física de su acción en el sistema. ¿Cuál es este efecto, y por qué se lo denota de esta forma?

Planteemos también el paso entre problemas de dominio acotado y los de no acotado. ¿Podemos aproximar las soluciones del sistema no acotado mediante las del acotado? Y de ser así, ¿qué condiciones de frontera en una dimensión funcionan mejor?

=Solución acotada vs Solución no acotada=
Demos por conocidos el uso de separación de variables para resolver problemas de autofunciones en caso acotado y el manejo de la solución fundamental junto con la convolución para el sistema no acotado. Para comparar todas las soluciones correspondientes, consideremos sistema <math> 1-</math>dimensional, con una gaussiana <math> e^{-x^2} </math> como condición inicial. Justificamos físicamente como el uso de un láser para calentar un punto de una tira fina de metal. Considerando que los fotones puedan desviarse como una distribución normal por el teorema central del límite, es un montaje experimental factible y regular. ¿Y qué condiciones de frontera?

Si consideramos de Dirichlet, a una tira de longitud <math> 2a </math> le colocaremos en los extremos un material que se mantenga siempre a la misma temperatura, pongamos que a una temperatura nula. El sistema resulta

<math>\quad
\begin{align}
\begin{cases} u_{t} - u_{xx} = 0 \quad &\forall t > 0, x\in [-a,a],\\
u(-a,t) = u(a,t) = 0 \quad &\forall t > 0, \\
u(x,0) = e^{-x^2} \quad &\forall x\in [-a,a].
\end{cases}
\end{align}
</math>

La solución al resolver el problema de autofunciones viene dada por

<math>\quad
u(x,t) = \sum_{n=1}^\infty A_{n} \cos\left(\frac{\pi}{a}(n - \frac{1}{2}) x\right) e^{-\dfrac{\pi^2}{a^2}\left(n - \frac{1}{2}\right)^2 t},
\quad \text{ donde } \quad
A_n = \frac{\int_{-a}^{a} u(x,0) \cos\left(\frac{\pi}{a}(n - \frac{1}{2}) x\right)dx}{\int_{-a}^{a} \cos^2\left(\frac{\pi}{a}(n - \frac{1}{2}) x\right)dx}.
</math>

Consideremos ahora condiciones de frontera Neumann, con un flujo nulo que físicamente viene de poner extremos aislantes. Para cierto <math> b > 0,</math> y no necesariamente igual que con <math> a </math>. Consideraremos la ecuación del calor

<math>\quad
\begin{align}
\begin{cases} u_{t} - u_{xx} = 0 \quad &\forall t > 0 , x\in [-b,b],\\
u_x(-b,t) = u_x(b,t) = 0 \quad &\forall t > 0, \\
u(x,0) = e^{-x^2} \quad &\forall x\in[-b,b].
\end{cases}
\end{align}
</math>

De nuevo podemos resolver el sistema, por ser homogéneo, mediante separación de variables, obteniendo

<math>\quad
u(x,t) = \sum_{m=0}^\infty B_{m} cos\left(\frac{m\pi}{b} x\right) e^{-\dfrac{m^2\pi^2}{b^2} t},\quad \text{ donde } \quad
B_m = \frac{\int_{-b}^{b} u(x,0) \cos\left(\frac{m\pi}{b} x\right)dx}{\int_{-b}^{b} \cos^2\left(\frac{m\pi}{b}x\right)dx}.
</math>

Un último detalle antes de proceder a los errores es la obtención simbólica de la solución no acotada.

<math>\quad
\begin{align}
\begin{cases} u_{t} - u_{xx} = 0 \quad &\forall t>0, x\in \mathbb{R},\\
u(x,0) = e^{-x^2} \quad &\forall x\in \mathbb{R}
\end{cases}
\end{align}
\quad
\Longrightarrow
\quad
u(x,t) = \frac{1}{\sqrt{4\pi t}}e^{\frac{-x^2}{4t}} * e^{-x^2} = \frac{e^{\frac{-x^2}{4t+1}}}{\sqrt{4t+1}}.
</math>

Queremos comprobar cuán bien aproximan las soluciones del caso acotado Dirichlet y Neumann a la obtenida por la solución fundamental. Para ello, escojamos un tiempo máximo <math>t_{max} </math> hasta el cuál tomaremos el error en norma <math> L^{\infty} </math>. Queremos estudiar:

* El efecto de aumentar el rango de definición con <math> a,b </math>.
* El efecto de aumentar el número de términos en las series de los sistemas acotados.

=== Análisis de los errores ===
[[Archivo:GrupoADMR Calor Errores1.jpg|750px|thumb|right| Representación gráfica de las soluciones del problema acotado (<math>a = b = 20</math>) y no acotado truncando las series de Fourier en su término <math>100</math>, para un tiempo máximo de <math>1</math>; y de sus errores. Nótese el orden de magnitud de los errores -pequeño-, los picos generados por los errores numéricos dados por la fórmula del trapecio y la suavidad en la frontera dada por las condiciones en esta.]]

Consideramos los parámetros, donde <math> n </math> es el término de la serie de Fourier donde truncamos:

*<math> a = b = 20, </math>
*<math> t_{max} = 1, </math>
*<math> n = 100. </math>

Con <math> t_{max} </math> fijo podemos observar que al aumentar <math>a, b</math> y <math>n</math> el error entre las soluciones acotadas y no acotada disminuye. Nótese que se estanca a partir de cierto <math>n</math>. De hecho, el error llega a aumentar a partir de ciertos <math> a,b </math> si truncamos demasiado pronto o demasiado tarde por la aproximación del trapecio para particiones del espacio demasiado bastas -no observado por pericia en código final-.

{|style="margin: 0 auto;"
| [[Archivo:GrupoADMR Calor Errores MaxD.jpeg|550px|thumb|upright|Con <math>t_{max}=1</math>, vamos obteniendo errores máximos en escala logarítmica con diferentes longitudes y truncamientos entre la solución del problema no acotado y el acotado Dirichlet. Los errores son bastante pequeños si se toman suficientes términos de la serie y una discretización espacial adecuada (en otro caso, aunque no esté representado en esta versión final, los errores divergen).]]
| [[Archivo:GrupoADMR Calor Errores MaxNe.jpeg|550px|thumb|upright|Análogo a la figura adyacente. Nótese que ambas son buenas aproximaciones con suficiente longitud y términos de la serie.]]
|}

Al ampliar el tiempo con la <math>n</math> y fijada <math>a=b=15</math>, al hacer más grande <math>t_{max}</math> los errores máximos no aumentan. En cambio, si tomamos el máximo error para diferentes tiempos vemos con suficiente precisión que no hay diferencia pero que para truncamientos precoces Neumann tiene ligeramente menor error.

{|style="margin: 0 auto;"
| [[Archivo:GrupoADMR Calor Errores DNo.jpeg|550px|thumb|upright|Con dominio fijo, el error en escala logarítmica entre solución no acotada y sistema Dirichlet decae incluso con el paso del tiempo. Esto encaja con que el estado estacionario de Dirichlet será la constante de temperatura nula, al igual que la convolución con la solución fundamental.]]
| [[Archivo:GrupoADMR Calor Errores NoNe.jpeg|550px|thumb|upright|Con condiciones análogas en frontera Neumann, el resultado es igualmente excelente para aproximar el caso no acotado. Aunque en el caso Neumann la solución estacionaria no sea nula, sino un reparto equitativo del calor inicial, cuando el rango es suficientemente amplio, la aproximación será adecuada.]]
|}

{|style="margin: 0 auto;"
| [[Archivo:GrupoADMR Calor Errores DNe.jpeg|550px|thumb|upright|Podemos ver además que ambos sistemas se asemejan cada vez más, por lo que ambas aproximaciones son igualmente útiles para el caso no acotado.]]
|}

=Coeficiente de difusión=

[[Archivo:Fundamental.gif|500px|thumb|right|Solución fundamental de la ecuación del calor unidimensional para distintos valores del coeficiente de difusión. La representación gráfica parte del tiempo inicial <math/> 0.01, <math/> esperando una delta de Dirac en el instante inicial. Nótese que en este tiempo, a mayor coeficiente de difusión más se ha difundido el calor en el sistema, lo que se nota especialmente en el valor máximo observado y su disminución cuanto más crece el coeficiente.]]

Interpretemos el valor de <math/> D <math/> aprovechando el montaje físico anterior. Numéricamente vemos que cuanto mayor es este coeficiente más rápido se difunde el calor, evidenciado en la caída de la solución en tiempos cada vez más cortos.

[[Archivo:Convoluciones_teorica.gif|500px|thumb|right|Soluciones de la ecuación del calor unidimensional para distintos valores del coeficiente de difusión.]]

Tanto las soluciones del caso unidimensional no acotado como las del acotado -que pueden reformularse considerando <math/> D <math/> en su exponencial, o aproximando por la sección anterior- podemos ver que un mayor coeficiente <math/> D <math/> acelera el decaimiento de la exponencial por lo que tendrá un efecto similar, justificándose su nombre.

= Códigos =

<source lang="matlab">
%%% Código generado junto con Chat GPT para la representación de la solución fundamental de la ecuación del calor unidimensional para distintos valores del coeficiente de difusión

clc
clear all
close all

% Vector de valores de D (coeficiente de difusión)
DD = [1,2,3,4,5,6,7,8,9,10];

% Nombre del archivo GIF de salida
output_gif = 'Fundamental.gif';

% Crear la figura
figura = figure(1);
grid on;
view(3);

% Bucle para iterar sobre cada valor de D
for i = 1:length(DD)
D = DD(i);

% Función fundamental de la ecuación del calor
phi_D = @(x, t) exp(-x.^2./(4*D*t))./sqrt(4*D*pi*t);

% Definición de la malla
xx = -1:0.01:1; % Rango de x
tt = 10^(-2):0.01:1; % Rango de t
[X, T] = meshgrid(xx, tt); % Malla de coordenadas
PHI_D = phi_D(X, T); % Evaluación de la función

% Graficar la superficie
clf;
surf(X, T, PHI_D);
shading flat
shading interp
colormap('jet')
xlabel('x');
ylabel('t');
zlabel('$\Phi_D(x,t)$', 'Interpreter', 'latex');
title(['Solución fundamental de la ecuación del calor, cuando D=', num2str(D)])
axis([-1, 1, 0, 1, 0, 3])
drawnow;

% Captura de la imagen para el GIF
frame = getframe(figura);
img = frame2im(frame);
[imind, cm] = rgb2ind(img, 256);

% Escribir la imagen en el archivo GIF
if i == 1
imwrite(imind, cm, output_gif, 'gif', 'LoopCount', Inf, 'DelayTime', 0.7);
else
imwrite(imind, cm, output_gif, 'gif', 'WriteMode', 'append', 'DelayTime', 0.7);
end

end
</source>

<source lang="matlab">
%%% Código generado junto con Chat GPT para la representación de la soluciones de la ecuación del calor unidimensional para distintos valores del coeficiente de difusión

clc
clear all
close all

% Vector de valores de D (coeficiente de difusión)
DD = [1,10,20,30,40,50,60,70,80,90,100];

% Función gaussiana
gaus = @(x) (exp(-x.^2));

% Nombre del archivo GIF de salida
output_gif = 'Convoluciones_teorica.gif';

% Crear la figura
figura = figure(1);
grid on;
view(3);

% Bucle para iterar sobre cada valor de D
for h = 1:length(DD)
D = DD(h);

% Funciones fundamentales y convolución
phi_D = @(x, t) exp(-x.^2./(4*D*t))./sqrt(4*D*pi*t);
convo = @(x,t) (exp(-x.^2./(4*D.*t+1))/sqrt(4*D.*t+1));

% Definición de la malla
xx = -10:0.1:10; % Rango de x
tt = 10^(-2):0.01:1; % Rango de t

% Inicialización de la matriz de convolución
conv = zeros(length(xx), length(tt));

% Cálculo de la convolución
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
val_conv(i,j) = convo(auxx, auxt);
end
end

% Graficar la superficie
clf;
hold on;
view(3);
surf(xx, tt, val_conv.');
xlabel('x');
ylabel('t');
zlabel('$\Phi_D(x,t)*N(x)$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap('jet');
title(['Soluciones en dominio no acotado, cuando D=', num2str(D)]);
axis([-10, 10, 0, 1, 0, 1]);
hold off;
drawnow;

% Captura de la imagen para el GIF
frame = getframe(figura);
img = frame2im(frame);
[imind, cm] = rgb2ind(img, 256);

% Escribir la imagen en el archivo GIF
if h == 1
imwrite(imind, cm, output_gif, 'gif', 'LoopCount', Inf, 'DelayTime', 0.7);
else
imwrite(imind, cm, output_gif, 'gif', 'WriteMode', 'append', 'DelayTime', 0.7);
end

end
</source>

<source lang="matlab">
%%% Código generado junto con Chat GPT para la representación de los errores entre las soluciones

function [Diff12,Difft1,Difft2] = errores2(n_max, a, b,t_max)

% Función como condición inicial: una función Gaussiana
gaus = @(x) (exp(-x.^2));
u0 = @(x) exp(-x.^2);
D = 1; % Coeficiente de difusión

% Funciones fundamentales y convolución para el modelo de difusión
phi_D = @(x, t) exp(-x.^2./(4*D*t))./sqrt(4*D*pi*t);
convo = @(x,t) (exp(-x.^2./(4*D.*t+1))/sqrt(4*D.*t+1));

% Definición de la malla (espacio y tiempo)
xx = -a:0.01:a;
tt = 10^(-2):0.01:t_max;
[X, T] = meshgrid(xx, tt);

% Inicialización de la matriz de convolución
conv = zeros(length(xx), length(tt));

% Cálculo de la convolución
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
val_conv(j,i) = convo(auxx, auxt);
end
end

% Inicialización de la matriz de resultados para la solución Dirichlet
res_dir = zeros(length(xx), length(tt));

% Cálculo de la solución directa
for n=1:n_max
A(n) = trapz(xx, u0(xx) .* cos((pi/a)*(n - 1/2) * xx)) /...
trapz(xx, cos((pi/a)*(n - 1/2) * xx).^2);
end

% Evaluación de la solución directa en el dominio espacio-temporal
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
u = 0;

% Sumar los términos de la serie de Fourier
for n=1:n_max
u = u + A(n) .* cos((pi/a)*(n - 1/2) * auxx) .*...
exp(-((pi^2 / a^2) * (n - 1/2)^2) .* auxt);
end
res_dir(i,j) = u;
end
end

% Inicialización de la matriz de resultados para la solución en frontera Neumann
res_neu = zeros(length(xx), length(tt));

% Cálculo de la solución en frontera Neumann
for n=1:n_max
B(n) = trapz(xx, u0(xx) .* cos((n*pi/b) * xx)) /...
trapz(xx, cos((n*pi/b) * xx).^2);
end

% Evaluación de la solución en frontera Neumann en el dominio espacio-temporal
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
u = trapz(xx, u0(xx))./(2.*b);

% Sumar los términos de la serie de Fourier
for n=1:n_max
u = u + B(n) .* cos((n*pi/b) * auxx) .* exp(-((n^2 * pi^2 / b^2)) .* auxt);
end
res_neu(i,j) = u;
end
end

% Transponer los resultados para compararlos con la convolución
U1 = res_dir.';
U2 = res_neu.';

% Calcular las diferencias entre las soluciones y la convolución
Diff12 = abs(U1 - U2);
Difft1 = abs(val_conv - U1);
Difft2 = abs(val_conv - U2);

end
</source>

<source lang="matlab">
%%% Código generado junto con Chat GPT para la representación de las soluciones de la ecuación del calor unidimensional junto con la comparación de sus errores

clc
clear all
close all

% Función como condicion inicial
% Se define la función inicial gausiana y la función u0 que se utilizarán más adelante para los cálculos.
gaus = @(x) (exp(-x.^2));
u0 = @(x) exp(-x.^2);
D = 1;

% Funciones fundamentales y convolución
% Se definen las funciones que describen la solución fundamental de la ecuación y la convolución utilizada.
phi_D = @(x, t) exp(-x.^2./(4*D*t))./sqrt(4*D*pi*t);
convo = @(x,t) (exp(-x.^2./(4*D.*t+1))/sqrt(4*D.*t+1));

% Definición de la malla
% Se definen los parámetros de la malla para las variables espaciales (xx) y temporales (tt), y se crea la malla 2D.
a = 100;
b = 100;
tmax = 1;
n_max = 1000;

xx = linspace(-a, a, 1000);
tt = linspace(10^(-2), tmax, 1000);
[X, T] = meshgrid(xx, tt);

% Inicialización de la matriz de convolución
% Se crea una matriz vacía para almacenar los valores de la convolución.
conv = zeros(length(xx), length(tt));

% Cálculo de la convolución
% Se calcula la convolución para cada punto en la malla de la variable temporal tt y espacial xx.
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
val_conv(i,j) = convo(auxx, auxt);
end
end

% Resolución de la ecuación con condiciones de frontera Dirichlet
% Se calculan los coeficientes A para la expansión en series de Fourier y se calcula la solución directa en cada paso temporal.
res_dir = zeros(length(xx), length(tt));

for n=1:n_max
A(n) = trapz(xx, u0(xx) .* cos((pi/a)*(n - 1/2) * xx)) /...
trapz(xx, cos((pi/a)*(n - 1/2) * xx).^2);
end

for j = 1:length(tt)
auxt = tt(j);

for i = 1:length(xx)
auxx = xx(i);
u = 0;

for n=1:n_max
u = u + A(n) .* cos((pi/a)*(n - 1/2) * auxx) .*...
exp(-((pi^2 / a^2) * (n - 1/2)^2) .* auxt);

end
res_dir(j,i) = u;
end
end

% Resolución de la ecuación con condiciones de frontera Neumann
% Se calculan los coeficientes B para la expansión en series de Fourier con las condiciones de frontera Neumann y se calcula la solución en cada paso temporal.
res_neu = zeros(length(xx), length(tt));

for n=1:n_max
B(n) = trapz(xx, u0(xx) .* cos((n*pi/b) * xx)) /...
trapz(xx, cos((n*pi/b) * xx).^2);
end

for j = 1:length(tt)
auxt = tt(j);

for i = 1:length(xx)
auxx = xx(i);
u = trapz(xx, u0(xx))./(2.*b);

for n=1:n_max
u = u + B(n) .* cos((n*pi/b) * auxx) .* exp(-((n^2 * pi^2 / b^2)) .* auxt);

end
res_neu(j,i) = u;
end
end

% Comparación entre las soluciones obtenidas
% Se calcula la diferencia entre las soluciones directa y Neumann, así como entre las soluciones y la convolución.
U1 = res_dir;
U2 = res_neu;

Diff12 = abs(U1 - U2);
Difft1 = abs(val_conv.' - U1);
Difft2 = abs(val_conv.' - U2);

% Visualización de los resultados
% Se muestran las superficies 3D de las soluciones y sus diferencias en gráficos para comparar los resultados obtenidos.
figure(1);

subplot(2,4,2);
view(3);
surf(X, T, U1);
xlabel('x');
ylabel('t');
zlabel('$u_1(x,t)$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca,'jet');
axis([-a,a,0,tmax,0,1])
title('$u_1(x,t)$', 'Interpreter', 'latex');

subplot(2,4,3);
view(3);
surf(X, T, U2);
xlabel('x');
ylabel('t');
zlabel('$u_2(x,t)$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca, 'jet');
axis([-b,b,0,tmax,0,1])
title('$u_2(x,t)$', 'Interpreter', 'latex');

subplot(2,4,4);
view(3);
surf(X, T, Diff12);
xlabel('x');
ylabel('t');
zlabel('$|u_1(x,t)-u_2(x,t)|$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca, 'hot');
title('$|u_1(x,t)-u_2(x,t)|$', 'Interpreter', 'latex');

subplot(2,4,5)
view(3);
surf(X, T, val_conv.');
xlabel('x');
ylabel('t');
zlabel('$\Phi_D(x,t)*N(x)$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca, 'jet');
axis([-a,a,0,tmax,0,1])
title('$\Phi_D(x,t)*N(x)$', 'Interpreter', 'latex');

subplot(2,4,6);
view(3);
surf(X, T, Difft1);
xlabel('x');
ylabel('t');
zlabel('$|\Phi_D(x,t)*N(x)-u_1(x,t)|$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca, 'hot');
title('$|\Phi_D(x,t)*N(x)-u_1(x,t)|$', 'Interpreter', 'latex');

subplot(2,4,7);
view(3);
surf(X, T, Difft2);
xlabel('x');
ylabel('t');
zlabel('$|\Phi_D(x,t)*N(x)-u_2(x,t)|$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca, 'hot');
title('$|\Phi_D(x,t)*N(x)-u_2(x,t)|$', 'Interpreter', 'latex');
</source>

= Referencias =

* Sandro Salsa, Gianmaria Verzini. Partial Differential Equations in Action From Modelling to Theory ([https://link.springer.com/book/10.1007/978-3-031-21853-8 [1<nowiki>]</nowiki>]).

[[Categoría:EDP]]
[[Categoría:EDP24/25]]

Ecuación del calor (ADMR)

2025-03-19T13:28:48Z

Daniel Rodríguez Calderón: /* Análisis de los errores */

{{ TrabajoED | Ecuación del calor (Grupo ADMR). | [[:Categoría:EDP|EDP]]|[[:Categoría:EDP24/25|2024-25]] |
*Ángel De Lucas Miranda,
*Daniel Rodríguez Calderón,
*Marcos Cabellos Hernández,
*Rafael Pascual Ortega.}}

=Introducción=

El desarrollo de la ecuación del calor se atribuye principalmente al matemático y físico francés Joseph Fourier, quien en 1822 publicó su obra ''Théorie analytique de la chaleur''. En este trabajo formuló una ecuación diferencial en derivadas parciales que describe la propagación del calor por difusión en un medio continuo. En esta ocasión, profundicemos en dos aspectos de este modelo.

Desde un punto de vista físico, estas ecuaciones permiten la modelización de fenómenos como el calentamiento de materiales o la difusión de solutos en fluidos estacionarios. Sin embargo, la influencia de la física es más profunda: de hecho, el coeficiente de difusión lleva su nombre por la interpretación física de su acción en el sistema. ¿Cuál es este efecto, y por qué se lo denota de esta forma?

Planteemos también el paso entre problemas de dominio acotado y los de no acotado. ¿Podemos aproximar las soluciones del sistema no acotado mediante las del acotado? Y de ser así, ¿qué condiciones de frontera en una dimensión funcionan mejor?

=Solución acotada vs Solución no acotada=
Demos por conocidos el uso de separación de variables para resolver problemas de autofunciones en caso acotado y el manejo de la solución fundamental junto con la convolución para el sistema no acotado. Para comparar todas las soluciones correspondientes, consideremos sistema <math> 1-</math>dimensional, con una gaussiana <math> e^{-x^2} </math> como condición inicial. Justificamos físicamente como el uso de un láser para calentar un punto de una tira fina de metal. Considerando que los fotones puedan desviarse como una distribución normal por el teorema central del límite, es un montaje experimental factible y regular. ¿Y qué condiciones de frontera?

Si consideramos de Dirichlet, a una tira de longitud <math> 2a </math> le colocaremos en los extremos un material que se mantenga siempre a la misma temperatura, pongamos que a una temperatura nula. El sistema resulta

<math>\quad
\begin{align}
\begin{cases} u_{t} - u_{xx} = 0 \quad &\forall t > 0, x\in [-a,a],\\
u(-a,t) = u(a,t) = 0 \quad &\forall t > 0, \\
u(x,0) = e^{-x^2} \quad &\forall x\in [-a,a].
\end{cases}
\end{align}
</math>

La solución al resolver el problema de autofunciones viene dada por

<math>\quad
u(x,t) = \sum_{n=1}^\infty A_{n} \cos\left(\frac{\pi}{a}(n - \frac{1}{2}) x\right) e^{-\dfrac{\pi^2}{a^2}\left(n - \frac{1}{2}\right)^2 t},
\quad \text{ donde } \quad
A_n = \frac{\int_{-a}^{a} u(x,0) \cos\left(\frac{\pi}{a}(n - \frac{1}{2}) x\right)dx}{\int_{-a}^{a} \cos^2\left(\frac{\pi}{a}(n - \frac{1}{2}) x\right)dx}.
</math>

Consideremos ahora condiciones de frontera Neumann, con un flujo nulo que físicamente viene de poner extremos aislantes. Para cierto <math> b > 0,</math> y no necesariamente igual que con <math> a </math>. Consideraremos la ecuación del calor

<math>\quad
\begin{align}
\begin{cases} u_{t} - u_{xx} = 0 \quad &\forall t > 0 , x\in [-b,b],\\
u_x(-b,t) = u_x(b,t) = 0 \quad &\forall t > 0, \\
u(x,0) = e^{-x^2} \quad &\forall x\in[-b,b].
\end{cases}
\end{align}
</math>

De nuevo podemos resolver el sistema, por ser homogéneo, mediante separación de variables, obteniendo

<math>\quad
u(x,t) = \sum_{m=0}^\infty B_{m} cos\left(\frac{m\pi}{b} x\right) e^{-\dfrac{m^2\pi^2}{b^2} t},\quad \text{ donde } \quad
B_m = \frac{\int_{-b}^{b} u(x,0) \cos\left(\frac{m\pi}{b} x\right)dx}{\int_{-b}^{b} \cos^2\left(\frac{m\pi}{b}x\right)dx}.
</math>

Un último detalle antes de proceder a los errores es la obtención simbólica de la solución no acotada.

<math>\quad
\begin{align}
\begin{cases} u_{t} - u_{xx} = 0 \quad &\forall t>0, x\in \mathbb{R},\\
u(x,0) = e^{-x^2} \quad &\forall x\in \mathbb{R}
\end{cases}
\end{align}
\quad
\Longrightarrow
\quad
u(x,t) = \frac{1}{\sqrt{4\pi t}}e^{\frac{-x^2}{4t}} * e^{-x^2} = \frac{e^{\frac{-x^2}{4t+1}}}{\sqrt{4t+1}}.
</math>

Queremos comprobar cuán bien aproximan las soluciones del caso acotado Dirichlet y Neumann a la obtenida por la solución fundamental. Para ello, escojamos un tiempo máximo <math>t_{max} </math> hasta el cuál tomaremos el error en norma <math> L^{\infty} </math>. Queremos estudiar:

* El efecto de aumentar el rango de definición con <math> a,b </math>.
* El efecto de aumentar el número de términos en las series de los sistemas acotados.

=== Análisis de los errores ===
[[Archivo:GrupoADMR Calor Errores1.jpg|750px|thumb|right| Representación gráfica de las soluciones del problema acotado (<math>a = b = 20</math>) y no acotado truncando las series de Fourier en su término <math>100</math>, para un tiempo máximo de <math>1</math>; y de sus errores. Nótese el orden de magnitud de los errores -pequeño-, los picos generados por los errores numéricos dados por la fórmula del trapecio y la suavidad en la frontera dada por las condiciones en esta.]]

Consideramos los parámetros, donde <math> n </math> es el término de la serie de Fourier donde truncamos:

*<math> a = b = 20, </math>
*<math> t_{max} = 1, </math>
*<math> n = 100. </math>

Con <math> t_{max} </math> fijo podemos observar que al aumentar <math>a, b</math> y <math>n</math> el error entre las soluciones acotadas y no acotada disminuye. Nótese que se estanca a partir de cierto <math>n</math>. De hecho, el error llega a aumentar a partir de ciertos <math> a,b </math> si truncamos demasiado pronto o demasiado tarde por la aproximación del trapecio para particiones del espacio demasiado bastas -no observado por pericia en código final-.

{|style="margin: 0 auto;"
| [[Archivo:GrupoADMR Calor Errores MaxD.jpeg|550px|thumb|upright|Con <math>t_{max}=1</math>, vamos obteniendo errores máximos en escala logarítmica con diferentes longitudes y truncamientos entre la solución del problema no acotado y el acotado Dirichlet. Los errores son bastante pequeños si se toman suficientes términos de la serie y una discretización espacial adecuada (en otro caso, aunque no esté representado en esta versión final, los errores divergen).]]
| [[Archivo:GrupoADMR Calor Errores MaxNe.jpeg|550px|thumb|upright|Análogo a la figura adyacente. Nótese que ambas son buenas aproximaciones con suficiente longitud y términos de la serie.]]
|}

Al ampliar el tiempo con la <math>n</math> y fijada <math>a=b=15</math>, al hacer más grande <math>t_{max}</math> los errores máximos no aumentan. En cambio, si tomamos el máximo error para diferentes tiempos vemos con suficiente precisión que no hay diferencia pero que para truncamientos precoces Neumann tiene ligeramente menor error.

{|style="margin: 0 auto;"
| [[Archivo:GrupoADMR Calor Errores DNo.jpeg|550px|thumb|upright|Con dominio fijo, el error en escala logarítmica entre solución no acotada y sistema Dirichlet decae incluso con el paso del tiempo. Esto encaja con que el estado estacionario de Dirichlet será la constante de temperatura nula, al igual que la convolución con la solución fundamental.]]
| [[Archivo:GrupoADMR Calor Errores NoNe.jpeg|550px|thumb|upright|Con condiciones análogas en frontera Neumann, el resultado es igualmente excelente para aproximar el caso no acotado. Aunque en el caso Neumann la solución estacionaria no sea nula, sino un reparto equitativo del calor inicial, cuando el rango es suficientemente amplio, la aproximación será adecuada.]]
|}

{|style="margin: 0 auto;"
| [[Archivo:GrupoADMR Calor Errores DNe.jpeg|550px|thumb|upright|Podemos ver además que ambos sistemas se asemejan cada vez más, por lo que ambas aproximaciones son igualmente útiles para el caso no acotado.]]
|}

=Coeficiente de difusión=

[[Archivo:Fundamental.gif|500px|thumb|right|Solución fundamental de la ecuación del calor unidimensional para distintos valores del coeficiente de difusión. La representación gráfica parte del tiempo inicial <math/> 0.01, <math/> esperando una delta de Dirac en el instante inicial. Nótese que en este tiempo, a mayor coeficiente de difusión más se ha difundido el calor en el sistema, lo que se nota especialmente en el valor máximo observado y su disminución cuanto más crece el coeficiente.]]

Interpretemos el valor de <math/> D <math/> aprovechando el montaje físico anterior. Numéricamente vemos que cuanto mayor es este coeficiente más rápido se difunde el calor, evidenciado en la caída de la solución en tiempos cada vez más cortos.

[[Archivo:Convoluciones_teorica.gif|500px|thumb|right|Soluciones de la ecuación del calor unidimensional para distintos valores del coeficiente de difusión.]]

En el caso unidimensional acotado con solución

<math>
u(x,t) = \sum_{n=1}^{\infty} A_n e^{-D (n\pi)^2 t} \sin(n\pi x).
</math>

Podemos ver que un mayor coeficiente <math/> D <math/> acelera el decaimiento de la exponencial por lo que tendrá un efecto similar justificando su nombre.

= Códigos =

<source lang="matlab">
%%% Código generado junto con Chat GPT para la representación de la solución fundamental de la ecuación del calor unidimensional para distintos valores del coeficiente de difusión

clc
clear all
close all

% Vector de valores de D (coeficiente de difusión)
DD = [1,2,3,4,5,6,7,8,9,10];

% Nombre del archivo GIF de salida
output_gif = 'Fundamental.gif';

% Crear la figura
figura = figure(1);
grid on;
view(3);

% Bucle para iterar sobre cada valor de D
for i = 1:length(DD)
D = DD(i);

% Función fundamental de la ecuación del calor
phi_D = @(x, t) exp(-x.^2./(4*D*t))./sqrt(4*D*pi*t);

% Definición de la malla
xx = -1:0.01:1; % Rango de x
tt = 10^(-2):0.01:1; % Rango de t
[X, T] = meshgrid(xx, tt); % Malla de coordenadas
PHI_D = phi_D(X, T); % Evaluación de la función

% Graficar la superficie
clf;
surf(X, T, PHI_D);
shading flat
shading interp
colormap('jet')
xlabel('x');
ylabel('t');
zlabel('$\Phi_D(x,t)$', 'Interpreter', 'latex');
title(['Solución fundamental de la ecuación del calor, cuando D=', num2str(D)])
axis([-1, 1, 0, 1, 0, 3])
drawnow;

% Captura de la imagen para el GIF
frame = getframe(figura);
img = frame2im(frame);
[imind, cm] = rgb2ind(img, 256);

% Escribir la imagen en el archivo GIF
if i == 1
imwrite(imind, cm, output_gif, 'gif', 'LoopCount', Inf, 'DelayTime', 0.7);
else
imwrite(imind, cm, output_gif, 'gif', 'WriteMode', 'append', 'DelayTime', 0.7);
end

end
</source>

<source lang="matlab">
%%% Código generado junto con Chat GPT para la representación de la soluciones de la ecuación del calor unidimensional para distintos valores del coeficiente de difusión

clc
clear all
close all

% Vector de valores de D (coeficiente de difusión)
DD = [1,10,20,30,40,50,60,70,80,90,100];

% Función gaussiana
gaus = @(x) (exp(-x.^2));

% Nombre del archivo GIF de salida
output_gif = 'Convoluciones_teorica.gif';

% Crear la figura
figura = figure(1);
grid on;
view(3);

% Bucle para iterar sobre cada valor de D
for h = 1:length(DD)
D = DD(h);

% Funciones fundamentales y convolución
phi_D = @(x, t) exp(-x.^2./(4*D*t))./sqrt(4*D*pi*t);
convo = @(x,t) (exp(-x.^2./(4*D.*t+1))/sqrt(4*D.*t+1));

% Definición de la malla
xx = -10:0.1:10; % Rango de x
tt = 10^(-2):0.01:1; % Rango de t

% Inicialización de la matriz de convolución
conv = zeros(length(xx), length(tt));

% Cálculo de la convolución
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
val_conv(i,j) = convo(auxx, auxt);
end
end

% Graficar la superficie
clf;
hold on;
view(3);
surf(xx, tt, val_conv.');
xlabel('x');
ylabel('t');
zlabel('$\Phi_D(x,t)*N(x)$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap('jet');
title(['Soluciones en dominio no acotado, cuando D=', num2str(D)]);
axis([-10, 10, 0, 1, 0, 1]);
hold off;
drawnow;

% Captura de la imagen para el GIF
frame = getframe(figura);
img = frame2im(frame);
[imind, cm] = rgb2ind(img, 256);

% Escribir la imagen en el archivo GIF
if h == 1
imwrite(imind, cm, output_gif, 'gif', 'LoopCount', Inf, 'DelayTime', 0.7);
else
imwrite(imind, cm, output_gif, 'gif', 'WriteMode', 'append', 'DelayTime', 0.7);
end

end
</source>

<source lang="matlab">
%%% Código generado junto con Chat GPT para la representación de los errores entre las soluciones

function [Diff12,Difft1,Difft2] = errores2(n_max, a, b,t_max)

% Función como condición inicial: una función Gaussiana
gaus = @(x) (exp(-x.^2));
u0 = @(x) exp(-x.^2);
D = 1; % Coeficiente de difusión

% Funciones fundamentales y convolución para el modelo de difusión
phi_D = @(x, t) exp(-x.^2./(4*D*t))./sqrt(4*D*pi*t);
convo = @(x,t) (exp(-x.^2./(4*D.*t+1))/sqrt(4*D.*t+1));

% Definición de la malla (espacio y tiempo)
xx = -a:0.01:a;
tt = 10^(-2):0.01:t_max;
[X, T] = meshgrid(xx, tt);

% Inicialización de la matriz de convolución
conv = zeros(length(xx), length(tt));

% Cálculo de la convolución
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
val_conv(j,i) = convo(auxx, auxt);
end
end

% Inicialización de la matriz de resultados para la solución Dirichlet
res_dir = zeros(length(xx), length(tt));

% Cálculo de la solución directa
for n=1:n_max
A(n) = trapz(xx, u0(xx) .* cos((pi/a)*(n - 1/2) * xx)) /...
trapz(xx, cos((pi/a)*(n - 1/2) * xx).^2);
end

% Evaluación de la solución directa en el dominio espacio-temporal
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
u = 0;

% Sumar los términos de la serie de Fourier
for n=1:n_max
u = u + A(n) .* cos((pi/a)*(n - 1/2) * auxx) .*...
exp(-((pi^2 / a^2) * (n - 1/2)^2) .* auxt);
end
res_dir(i,j) = u;
end
end

% Inicialización de la matriz de resultados para la solución en frontera Neumann
res_neu = zeros(length(xx), length(tt));

% Cálculo de la solución en frontera Neumann
for n=1:n_max
B(n) = trapz(xx, u0(xx) .* cos((n*pi/b) * xx)) /...
trapz(xx, cos((n*pi/b) * xx).^2);
end

% Evaluación de la solución en frontera Neumann en el dominio espacio-temporal
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
u = trapz(xx, u0(xx))./(2.*b);

% Sumar los términos de la serie de Fourier
for n=1:n_max
u = u + B(n) .* cos((n*pi/b) * auxx) .* exp(-((n^2 * pi^2 / b^2)) .* auxt);
end
res_neu(i,j) = u;
end
end

% Transponer los resultados para compararlos con la convolución
U1 = res_dir.';
U2 = res_neu.';

% Calcular las diferencias entre las soluciones y la convolución
Diff12 = abs(U1 - U2);
Difft1 = abs(val_conv - U1);
Difft2 = abs(val_conv - U2);

end
</source>

<source lang="matlab">
%%% Código generado junto con Chat GPT para la representación de las soluciones de la ecuación del calor unidimensional junto con la comparación de sus errores

clc
clear all
close all

% Función como condicion inicial
% Se define la función inicial gausiana y la función u0 que se utilizarán más adelante para los cálculos.
gaus = @(x) (exp(-x.^2));
u0 = @(x) exp(-x.^2);
D = 1;

% Funciones fundamentales y convolución
% Se definen las funciones que describen la solución fundamental de la ecuación y la convolución utilizada.
phi_D = @(x, t) exp(-x.^2./(4*D*t))./sqrt(4*D*pi*t);
convo = @(x,t) (exp(-x.^2./(4*D.*t+1))/sqrt(4*D.*t+1));

% Definición de la malla
% Se definen los parámetros de la malla para las variables espaciales (xx) y temporales (tt), y se crea la malla 2D.
a = 100;
b = 100;
tmax = 1;
n_max = 1000;

xx = linspace(-a, a, 1000);
tt = linspace(10^(-2), tmax, 1000);
[X, T] = meshgrid(xx, tt);

% Inicialización de la matriz de convolución
% Se crea una matriz vacía para almacenar los valores de la convolución.
conv = zeros(length(xx), length(tt));

% Cálculo de la convolución
% Se calcula la convolución para cada punto en la malla de la variable temporal tt y espacial xx.
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
val_conv(i,j) = convo(auxx, auxt);
end
end

% Resolución de la ecuación con condiciones de frontera Dirichlet
% Se calculan los coeficientes A para la expansión en series de Fourier y se calcula la solución directa en cada paso temporal.
res_dir = zeros(length(xx), length(tt));

for n=1:n_max
A(n) = trapz(xx, u0(xx) .* cos((pi/a)*(n - 1/2) * xx)) /...
trapz(xx, cos((pi/a)*(n - 1/2) * xx).^2);
end

for j = 1:length(tt)
auxt = tt(j);

for i = 1:length(xx)
auxx = xx(i);
u = 0;

for n=1:n_max
u = u + A(n) .* cos((pi/a)*(n - 1/2) * auxx) .*...
exp(-((pi^2 / a^2) * (n - 1/2)^2) .* auxt);

end
res_dir(j,i) = u;
end
end

% Resolución de la ecuación con condiciones de frontera Neumann
% Se calculan los coeficientes B para la expansión en series de Fourier con las condiciones de frontera Neumann y se calcula la solución en cada paso temporal.
res_neu = zeros(length(xx), length(tt));

for n=1:n_max
B(n) = trapz(xx, u0(xx) .* cos((n*pi/b) * xx)) /...
trapz(xx, cos((n*pi/b) * xx).^2);
end

for j = 1:length(tt)
auxt = tt(j);

for i = 1:length(xx)
auxx = xx(i);
u = trapz(xx, u0(xx))./(2.*b);

for n=1:n_max
u = u + B(n) .* cos((n*pi/b) * auxx) .* exp(-((n^2 * pi^2 / b^2)) .* auxt);

end
res_neu(j,i) = u;
end
end

% Comparación entre las soluciones obtenidas
% Se calcula la diferencia entre las soluciones directa y Neumann, así como entre las soluciones y la convolución.
U1 = res_dir;
U2 = res_neu;

Diff12 = abs(U1 - U2);
Difft1 = abs(val_conv.' - U1);
Difft2 = abs(val_conv.' - U2);

% Visualización de los resultados
% Se muestran las superficies 3D de las soluciones y sus diferencias en gráficos para comparar los resultados obtenidos.
figure(1);

subplot(2,4,2);
view(3);
surf(X, T, U1);
xlabel('x');
ylabel('t');
zlabel('$u_1(x,t)$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca,'jet');
axis([-a,a,0,tmax,0,1])
title('$u_1(x,t)$', 'Interpreter', 'latex');

subplot(2,4,3);
view(3);
surf(X, T, U2);
xlabel('x');
ylabel('t');
zlabel('$u_2(x,t)$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca, 'jet');
axis([-b,b,0,tmax,0,1])
title('$u_2(x,t)$', 'Interpreter', 'latex');

subplot(2,4,4);
view(3);
surf(X, T, Diff12);
xlabel('x');
ylabel('t');
zlabel('$|u_1(x,t)-u_2(x,t)|$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca, 'hot');
title('$|u_1(x,t)-u_2(x,t)|$', 'Interpreter', 'latex');

subplot(2,4,5)
view(3);
surf(X, T, val_conv.');
xlabel('x');
ylabel('t');
zlabel('$\Phi_D(x,t)*N(x)$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca, 'jet');
axis([-a,a,0,tmax,0,1])
title('$\Phi_D(x,t)*N(x)$', 'Interpreter', 'latex');

subplot(2,4,6);
view(3);
surf(X, T, Difft1);
xlabel('x');
ylabel('t');
zlabel('$|\Phi_D(x,t)*N(x)-u_1(x,t)|$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca, 'hot');
title('$|\Phi_D(x,t)*N(x)-u_1(x,t)|$', 'Interpreter', 'latex');

subplot(2,4,7);
view(3);
surf(X, T, Difft2);
xlabel('x');
ylabel('t');
zlabel('$|\Phi_D(x,t)*N(x)-u_2(x,t)|$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca, 'hot');
title('$|\Phi_D(x,t)*N(x)-u_2(x,t)|$', 'Interpreter', 'latex');
</source>

= Referencias =

* Sandro Salsa, Gianmaria Verzini. Partial Differential Equations in Action From Modelling to Theory ([https://link.springer.com/book/10.1007/978-3-031-21853-8 [1<nowiki>]</nowiki>]).

[[Categoría:EDP]]
[[Categoría:EDP24/25]]

Ecuación del calor (ADMR)

2025-03-19T13:25:19Z

Daniel Rodríguez Calderón: /* Coeficiente de difusión */

{{ TrabajoED | Ecuación del calor (Grupo ADMR). | [[:Categoría:EDP|EDP]]|[[:Categoría:EDP24/25|2024-25]] |
*Ángel De Lucas Miranda,
*Daniel Rodríguez Calderón,
*Marcos Cabellos Hernández,
*Rafael Pascual Ortega.}}

=Introducción=

El desarrollo de la ecuación del calor se atribuye principalmente al matemático y físico francés Joseph Fourier, quien en 1822 publicó su obra ''Théorie analytique de la chaleur''. En este trabajo formuló una ecuación diferencial en derivadas parciales que describe la propagación del calor por difusión en un medio continuo. En esta ocasión, profundicemos en dos aspectos de este modelo.

Desde un punto de vista físico, estas ecuaciones permiten la modelización de fenómenos como el calentamiento de materiales o la difusión de solutos en fluidos estacionarios. Sin embargo, la influencia de la física es más profunda: de hecho, el coeficiente de difusión lleva su nombre por la interpretación física de su acción en el sistema. ¿Cuál es este efecto, y por qué se lo denota de esta forma?

Planteemos también el paso entre problemas de dominio acotado y los de no acotado. ¿Podemos aproximar las soluciones del sistema no acotado mediante las del acotado? Y de ser así, ¿qué condiciones de frontera en una dimensión funcionan mejor?

=Solución acotada vs Solución no acotada=
Demos por conocidos el uso de separación de variables para resolver problemas de autofunciones en caso acotado y el manejo de la solución fundamental junto con la convolución para el sistema no acotado. Para comparar todas las soluciones correspondientes, consideremos sistema <math> 1-</math>dimensional, con una gaussiana <math> e^{-x^2} </math> como condición inicial. Justificamos físicamente como el uso de un láser para calentar un punto de una tira fina de metal. Considerando que los fotones puedan desviarse como una distribución normal por el teorema central del límite, es un montaje experimental factible y regular. ¿Y qué condiciones de frontera?

Si consideramos de Dirichlet, a una tira de longitud <math> 2a </math> le colocaremos en los extremos un material que se mantenga siempre a la misma temperatura, pongamos que a una temperatura nula. El sistema resulta

<math>\quad
\begin{align}
\begin{cases} u_{t} - u_{xx} = 0 \quad &\forall t > 0, x\in [-a,a],\\
u(-a,t) = u(a,t) = 0 \quad &\forall t > 0, \\
u(x,0) = e^{-x^2} \quad &\forall x\in [-a,a].
\end{cases}
\end{align}
</math>

La solución al resolver el problema de autofunciones viene dada por

<math>\quad
u(x,t) = \sum_{n=1}^\infty A_{n} \cos\left(\frac{\pi}{a}(n - \frac{1}{2}) x\right) e^{-\dfrac{\pi^2}{a^2}\left(n - \frac{1}{2}\right)^2 t},
\quad \text{ donde } \quad
A_n = \frac{\int_{-a}^{a} u(x,0) \cos\left(\frac{\pi}{a}(n - \frac{1}{2}) x\right)dx}{\int_{-a}^{a} \cos^2\left(\frac{\pi}{a}(n - \frac{1}{2}) x\right)dx}.
</math>

Consideremos ahora condiciones de frontera Neumann, con un flujo nulo que físicamente viene de poner extremos aislantes. Para cierto <math> b > 0,</math> y no necesariamente igual que con <math> a </math>. Consideraremos la ecuación del calor

<math>\quad
\begin{align}
\begin{cases} u_{t} - u_{xx} = 0 \quad &\forall t > 0 , x\in [-b,b],\\
u_x(-b,t) = u_x(b,t) = 0 \quad &\forall t > 0, \\
u(x,0) = e^{-x^2} \quad &\forall x\in[-b,b].
\end{cases}
\end{align}
</math>

De nuevo podemos resolver el sistema, por ser homogéneo, mediante separación de variables, obteniendo

<math>\quad
u(x,t) = \sum_{m=0}^\infty B_{m} cos\left(\frac{m\pi}{b} x\right) e^{-\dfrac{m^2\pi^2}{b^2} t},\quad \text{ donde } \quad
B_m = \frac{\int_{-b}^{b} u(x,0) \cos\left(\frac{m\pi}{b} x\right)dx}{\int_{-b}^{b} \cos^2\left(\frac{m\pi}{b}x\right)dx}.
</math>

Un último detalle antes de proceder a los errores es la obtención simbólica de la solución no acotada.

<math>\quad
\begin{align}
\begin{cases} u_{t} - u_{xx} = 0 \quad &\forall t>0, x\in \mathbb{R},\\
u(x,0) = e^{-x^2} \quad &\forall x\in \mathbb{R}
\end{cases}
\end{align}
\quad
\Longrightarrow
\quad
u(x,t) = \frac{1}{\sqrt{4\pi t}}e^{\frac{-x^2}{4t}} * e^{-x^2} = \frac{e^{\frac{-x^2}{4t+1}}}{\sqrt{4t+1}}.
</math>

Queremos comprobar cuán bien aproximan las soluciones del caso acotado Dirichlet y Neumann a la obtenida por la solución fundamental. Para ello, escojamos un tiempo máximo <math>t_{max} </math> hasta el cuál tomaremos el error en norma <math> L^{\infty} </math>. Queremos estudiar:

* El efecto de aumentar el rango de definición con <math> a,b </math>.
* El efecto de aumentar el número de términos en las series de los sistemas acotados.

=== Análisis de los errores ===
[[Archivo:GrupoADMR Calor Errores1.jpg|550px|thumb|right| Representación gráfica de las soluciones del problema acotado (<math>a = b = 20</math>) y no acotado truncando las series de Fourier en su término <math>100</math>, para un tiempo máximo de <math>1</math>; y de sus errores. Nótese el orden de magnitud de los errores -pequeño-, los picos generados por los errores numéricos dados por la fórmula del trapecio y la suavidad en la frontera dada por las condiciones en esta.]]

Consideramos los parámetros, donde <math> n </math> es el término de la serie de Fourier donde truncamos:

*<math> a = b = 20, </math>
*<math> t_{max} = 1, </math>
*<math> n = 100. </math>

Con <math> t_{max} </math> fijo podemos observar que al aumentar <math>a, b</math> y <math>n</math> el error entre las soluciones acotadas y no acotada disminuye. Nótese que se estanca a partir de cierto <math>n</math>. De hecho, el error llega a aumentar a partir de ciertos <math> a,b </math> si truncamos demasiado pronto o demasiado tarde por la aproximación del trapecio para particiones del espacio demasiado bastas -no observado por pericia en código final-.

{|style="margin: 0 auto;"
| [[Archivo:GrupoADMR Calor Errores MaxD.jpeg|250px|thumb|upright|Con <math>t_{max}=1</math>, vamos obteniendo errores máximos en escala logarítmica con diferentes longitudes y truncamientos entre la solución del problema no acotado y el acotado Dirichlet. Los errores son bastante pequeños si se toman suficientes términos de la serie y una discretización espacial adecuada (en otro caso, aunque no esté representado en esta versión final, los errores divergen)]]
| [[Archivo:GrupoADMR Calor Errores MaxNe.jpeg|250px|thumb|upright|Análogo a la figura adyacente. Nótese que ambas son buenas aproximaciones con suficiente longitud y términos de la serie.]]
|}

Al ampliar el tiempo con la <math>n</math> y fijada <math>a=b=15</math>, al hacer más grande <math>t_{max}</math> los errores máximos no aumentan. En cambio, si tomamos el máximo error para diferentes tiempos vemos con suficiente precisión que no hay diferencia pero que para truncamientos precoces Neumann tiene ligeramente menor error.

{|style="margin: 0 auto;"
| [[Archivo:GrupoADMR Calor Errores DNo.jpeg|250px|thumb|upright|Con dominio fijo, el error en escala logarítmica entre solución no acotada y sistema Dirichlet decae incluso con el paso del tiempo. Esto encaja con que el estado estacionario de Dirichlet será la constante de temperatura nula, al igual que la convolución con la solución fundamental.]]
| [[Archivo:GrupoADMR Calor Errores NoNe.jpeg|250px|thumb|upright|Con condiciones análogas en frontera Neumann, el resultado es igualmente excelente para aproximar el caso no acotado. Aunque en el caso Neumann la solución estacionaria no sea nula, sino un reparto equitativo del calor inicial, cuando el rango es suficientemente amplio, la aproximación será adecuada.]]
| [[Archivo:GrupoADMR Calor Errores DNe.jpeg|250px|thumb|upright|Podemos ver además que ambos sistemas se asemejan cada vez más, por lo que ambas aproximaciones son igualmente útiles para el caso no acotado.]]
|}

=Coeficiente de difusión=

[[Archivo:Fundamental.gif|500px|thumb|right|Solución fundamental de la ecuación del calor unidimensional para distintos valores del coeficiente de difusión. La representación gráfica parte del tiempo inicial <math/> 0.01, <math/> esperando una delta de Dirac en el instante inicial. Nótese que en este tiempo, a mayor coeficiente de difusión más se ha difundido el calor en el sistema, lo que se nota especialmente en el valor máximo observado y su disminución cuanto más crece el coeficiente.]]

Interpretemos el valor de <math/> D <math/> aprovechando el montaje físico anterior. Numéricamente vemos que cuanto mayor es este coeficiente más rápido se difunde el calor, evidenciado en la caída de la solución en tiempos cada vez más cortos.

[[Archivo:Convoluciones_teorica.gif|500px|thumb|right|Soluciones de la ecuación del calor unidimensional para distintos valores del coeficiente de difusión.]]

En el caso unidimensional acotado con solución

<math>
u(x,t) = \sum_{n=1}^{\infty} A_n e^{-D (n\pi)^2 t} \sin(n\pi x).
</math>

Podemos ver que un mayor coeficiente <math/> D <math/> acelera el decaimiento de la exponencial por lo que tendrá un efecto similar justificando su nombre.

= Códigos =

<source lang="matlab">
%%% Código generado junto con Chat GPT para la representación de la solución fundamental de la ecuación del calor unidimensional para distintos valores del coeficiente de difusión

clc
clear all
close all

% Vector de valores de D (coeficiente de difusión)
DD = [1,2,3,4,5,6,7,8,9,10];

% Nombre del archivo GIF de salida
output_gif = 'Fundamental.gif';

% Crear la figura
figura = figure(1);
grid on;
view(3);

% Bucle para iterar sobre cada valor de D
for i = 1:length(DD)
D = DD(i);

% Función fundamental de la ecuación del calor
phi_D = @(x, t) exp(-x.^2./(4*D*t))./sqrt(4*D*pi*t);

% Definición de la malla
xx = -1:0.01:1; % Rango de x
tt = 10^(-2):0.01:1; % Rango de t
[X, T] = meshgrid(xx, tt); % Malla de coordenadas
PHI_D = phi_D(X, T); % Evaluación de la función

% Graficar la superficie
clf;
surf(X, T, PHI_D);
shading flat
shading interp
colormap('jet')
xlabel('x');
ylabel('t');
zlabel('$\Phi_D(x,t)$', 'Interpreter', 'latex');
title(['Solución fundamental de la ecuación del calor, cuando D=', num2str(D)])
axis([-1, 1, 0, 1, 0, 3])
drawnow;

% Captura de la imagen para el GIF
frame = getframe(figura);
img = frame2im(frame);
[imind, cm] = rgb2ind(img, 256);

% Escribir la imagen en el archivo GIF
if i == 1
imwrite(imind, cm, output_gif, 'gif', 'LoopCount', Inf, 'DelayTime', 0.7);
else
imwrite(imind, cm, output_gif, 'gif', 'WriteMode', 'append', 'DelayTime', 0.7);
end

end
</source>

<source lang="matlab">
%%% Código generado junto con Chat GPT para la representación de la soluciones de la ecuación del calor unidimensional para distintos valores del coeficiente de difusión

clc
clear all
close all

% Vector de valores de D (coeficiente de difusión)
DD = [1,10,20,30,40,50,60,70,80,90,100];

% Función gaussiana
gaus = @(x) (exp(-x.^2));

% Nombre del archivo GIF de salida
output_gif = 'Convoluciones_teorica.gif';

% Crear la figura
figura = figure(1);
grid on;
view(3);

% Bucle para iterar sobre cada valor de D
for h = 1:length(DD)
D = DD(h);

% Funciones fundamentales y convolución
phi_D = @(x, t) exp(-x.^2./(4*D*t))./sqrt(4*D*pi*t);
convo = @(x,t) (exp(-x.^2./(4*D.*t+1))/sqrt(4*D.*t+1));

% Definición de la malla
xx = -10:0.1:10; % Rango de x
tt = 10^(-2):0.01:1; % Rango de t

% Inicialización de la matriz de convolución
conv = zeros(length(xx), length(tt));

% Cálculo de la convolución
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
val_conv(i,j) = convo(auxx, auxt);
end
end

% Graficar la superficie
clf;
hold on;
view(3);
surf(xx, tt, val_conv.');
xlabel('x');
ylabel('t');
zlabel('$\Phi_D(x,t)*N(x)$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap('jet');
title(['Soluciones en dominio no acotado, cuando D=', num2str(D)]);
axis([-10, 10, 0, 1, 0, 1]);
hold off;
drawnow;

% Captura de la imagen para el GIF
frame = getframe(figura);
img = frame2im(frame);
[imind, cm] = rgb2ind(img, 256);

% Escribir la imagen en el archivo GIF
if h == 1
imwrite(imind, cm, output_gif, 'gif', 'LoopCount', Inf, 'DelayTime', 0.7);
else
imwrite(imind, cm, output_gif, 'gif', 'WriteMode', 'append', 'DelayTime', 0.7);
end

end
</source>

<source lang="matlab">
%%% Código generado junto con Chat GPT para la representación de los errores entre las soluciones

function [Diff12,Difft1,Difft2] = errores2(n_max, a, b,t_max)

% Función como condición inicial: una función Gaussiana
gaus = @(x) (exp(-x.^2));
u0 = @(x) exp(-x.^2);
D = 1; % Coeficiente de difusión

% Funciones fundamentales y convolución para el modelo de difusión
phi_D = @(x, t) exp(-x.^2./(4*D*t))./sqrt(4*D*pi*t);
convo = @(x,t) (exp(-x.^2./(4*D.*t+1))/sqrt(4*D.*t+1));

% Definición de la malla (espacio y tiempo)
xx = -a:0.01:a;
tt = 10^(-2):0.01:t_max;
[X, T] = meshgrid(xx, tt);

% Inicialización de la matriz de convolución
conv = zeros(length(xx), length(tt));

% Cálculo de la convolución
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
val_conv(j,i) = convo(auxx, auxt);
end
end

% Inicialización de la matriz de resultados para la solución Dirichlet
res_dir = zeros(length(xx), length(tt));

% Cálculo de la solución directa
for n=1:n_max
A(n) = trapz(xx, u0(xx) .* cos((pi/a)*(n - 1/2) * xx)) /...
trapz(xx, cos((pi/a)*(n - 1/2) * xx).^2);
end

% Evaluación de la solución directa en el dominio espacio-temporal
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
u = 0;

% Sumar los términos de la serie de Fourier
for n=1:n_max
u = u + A(n) .* cos((pi/a)*(n - 1/2) * auxx) .*...
exp(-((pi^2 / a^2) * (n - 1/2)^2) .* auxt);
end
res_dir(i,j) = u;
end
end

% Inicialización de la matriz de resultados para la solución en frontera Neumann
res_neu = zeros(length(xx), length(tt));

% Cálculo de la solución en frontera Neumann
for n=1:n_max
B(n) = trapz(xx, u0(xx) .* cos((n*pi/b) * xx)) /...
trapz(xx, cos((n*pi/b) * xx).^2);
end

% Evaluación de la solución en frontera Neumann en el dominio espacio-temporal
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
u = trapz(xx, u0(xx))./(2.*b);

% Sumar los términos de la serie de Fourier
for n=1:n_max
u = u + B(n) .* cos((n*pi/b) * auxx) .* exp(-((n^2 * pi^2 / b^2)) .* auxt);
end
res_neu(i,j) = u;
end
end

% Transponer los resultados para compararlos con la convolución
U1 = res_dir.';
U2 = res_neu.';

% Calcular las diferencias entre las soluciones y la convolución
Diff12 = abs(U1 - U2);
Difft1 = abs(val_conv - U1);
Difft2 = abs(val_conv - U2);

end
</source>

<source lang="matlab">
%%% Código generado junto con Chat GPT para la representación de las soluciones de la ecuación del calor unidimensional junto con la comparación de sus errores

clc
clear all
close all

% Función como condicion inicial
% Se define la función inicial gausiana y la función u0 que se utilizarán más adelante para los cálculos.
gaus = @(x) (exp(-x.^2));
u0 = @(x) exp(-x.^2);
D = 1;

% Funciones fundamentales y convolución
% Se definen las funciones que describen la solución fundamental de la ecuación y la convolución utilizada.
phi_D = @(x, t) exp(-x.^2./(4*D*t))./sqrt(4*D*pi*t);
convo = @(x,t) (exp(-x.^2./(4*D.*t+1))/sqrt(4*D.*t+1));

% Definición de la malla
% Se definen los parámetros de la malla para las variables espaciales (xx) y temporales (tt), y se crea la malla 2D.
a = 100;
b = 100;
tmax = 1;
n_max = 1000;

xx = linspace(-a, a, 1000);
tt = linspace(10^(-2), tmax, 1000);
[X, T] = meshgrid(xx, tt);

% Inicialización de la matriz de convolución
% Se crea una matriz vacía para almacenar los valores de la convolución.
conv = zeros(length(xx), length(tt));

% Cálculo de la convolución
% Se calcula la convolución para cada punto en la malla de la variable temporal tt y espacial xx.
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
val_conv(i,j) = convo(auxx, auxt);
end
end

% Resolución de la ecuación con condiciones de frontera Dirichlet
% Se calculan los coeficientes A para la expansión en series de Fourier y se calcula la solución directa en cada paso temporal.
res_dir = zeros(length(xx), length(tt));

for n=1:n_max
A(n) = trapz(xx, u0(xx) .* cos((pi/a)*(n - 1/2) * xx)) /...
trapz(xx, cos((pi/a)*(n - 1/2) * xx).^2);
end

for j = 1:length(tt)
auxt = tt(j);

for i = 1:length(xx)
auxx = xx(i);
u = 0;

for n=1:n_max
u = u + A(n) .* cos((pi/a)*(n - 1/2) * auxx) .*...
exp(-((pi^2 / a^2) * (n - 1/2)^2) .* auxt);

end
res_dir(j,i) = u;
end
end

% Resolución de la ecuación con condiciones de frontera Neumann
% Se calculan los coeficientes B para la expansión en series de Fourier con las condiciones de frontera Neumann y se calcula la solución en cada paso temporal.
res_neu = zeros(length(xx), length(tt));

for n=1:n_max
B(n) = trapz(xx, u0(xx) .* cos((n*pi/b) * xx)) /...
trapz(xx, cos((n*pi/b) * xx).^2);
end

for j = 1:length(tt)
auxt = tt(j);

for i = 1:length(xx)
auxx = xx(i);
u = trapz(xx, u0(xx))./(2.*b);

for n=1:n_max
u = u + B(n) .* cos((n*pi/b) * auxx) .* exp(-((n^2 * pi^2 / b^2)) .* auxt);

end
res_neu(j,i) = u;
end
end

% Comparación entre las soluciones obtenidas
% Se calcula la diferencia entre las soluciones directa y Neumann, así como entre las soluciones y la convolución.
U1 = res_dir;
U2 = res_neu;

Diff12 = abs(U1 - U2);
Difft1 = abs(val_conv.' - U1);
Difft2 = abs(val_conv.' - U2);

% Visualización de los resultados
% Se muestran las superficies 3D de las soluciones y sus diferencias en gráficos para comparar los resultados obtenidos.
figure(1);

subplot(2,4,2);
view(3);
surf(X, T, U1);
xlabel('x');
ylabel('t');
zlabel('$u_1(x,t)$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca,'jet');
axis([-a,a,0,tmax,0,1])
title('$u_1(x,t)$', 'Interpreter', 'latex');

subplot(2,4,3);
view(3);
surf(X, T, U2);
xlabel('x');
ylabel('t');
zlabel('$u_2(x,t)$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca, 'jet');
axis([-b,b,0,tmax,0,1])
title('$u_2(x,t)$', 'Interpreter', 'latex');

subplot(2,4,4);
view(3);
surf(X, T, Diff12);
xlabel('x');
ylabel('t');
zlabel('$|u_1(x,t)-u_2(x,t)|$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca, 'hot');
title('$|u_1(x,t)-u_2(x,t)|$', 'Interpreter', 'latex');

subplot(2,4,5)
view(3);
surf(X, T, val_conv.');
xlabel('x');
ylabel('t');
zlabel('$\Phi_D(x,t)*N(x)$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca, 'jet');
axis([-a,a,0,tmax,0,1])
title('$\Phi_D(x,t)*N(x)$', 'Interpreter', 'latex');

subplot(2,4,6);
view(3);
surf(X, T, Difft1);
xlabel('x');
ylabel('t');
zlabel('$|\Phi_D(x,t)*N(x)-u_1(x,t)|$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca, 'hot');
title('$|\Phi_D(x,t)*N(x)-u_1(x,t)|$', 'Interpreter', 'latex');

subplot(2,4,7);
view(3);
surf(X, T, Difft2);
xlabel('x');
ylabel('t');
zlabel('$|\Phi_D(x,t)*N(x)-u_2(x,t)|$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca, 'hot');
title('$|\Phi_D(x,t)*N(x)-u_2(x,t)|$', 'Interpreter', 'latex');
</source>

= Referencias =

* Sandro Salsa, Gianmaria Verzini. Partial Differential Equations in Action From Modelling to Theory ([https://link.springer.com/book/10.1007/978-3-031-21853-8 [1<nowiki>]</nowiki>]).

[[Categoría:EDP]]
[[Categoría:EDP24/25]]

Ecuación del calor (ADMR)

2025-03-19T07:38:58Z

Daniel Rodríguez Calderón: /* Códigos */

{{ TrabajoED | Ecuación del calor (Grupo ADMR). | [[:Categoría:EDP|EDP]]|[[:Categoría:EDP24/25|2024-25]] |
*Ángel De Lucas Miranda,
*Daniel Rodríguez Calderón,
*Marcos Cabellos Hernández,
*Rafael Pascual Ortega.}}

=Introducción=

El desarrollo de la ecuación del calor se atribuye principalmente al matemático y físico francés Joseph Fourier, quien en 1822 publicó su obra ''Théorie analytique de la chaleur''. En este trabajo formuló una ecuación diferencial en derivadas parciales que describe la propagación del calor por difusión en un medio continuo. En esta ocasión, profundicemos en dos aspectos de este modelo.

Desde un punto de vista físico, estas ecuaciones permiten la modelización de fenómenos como el calentamiento de materiales o la difusión de solutos en fluidos estacionarios. Sin embargo, la influencia de la física es más profunda: de hecho, el coeficiente de difusión lleva su nombre por la interpretación física de su acción en el sistema. ¿Cuál es este efecto, y por qué se lo denota de esta forma?

Planteemos también el paso entre problemas de dominio acotado y los de no acotado. ¿Podemos aproximar las soluciones del sistema no acotado mediante las del acotado? Y de ser así, ¿qué condiciones de frontera en una dimensión funcionan mejor?

=Solución acotada vs Solución no acotada=
Demos por conocidos el uso de separación de variables para resolver problemas de autofunciones en caso acotado y el manejo de la solución fundamental junto con la convolución para el sistema no acotado. Para comparar todas las soluciones correspondientes, consideremos sistema <math> 1-</math>dimensional, con una gaussiana <math> e^{-x^2} </math> como condición inicial. Justificamos físicamente como el uso de un láser para calentar un punto de una tira fina de metal. Considerando que los fotones puedan desviarse como una distribución normal por el teorema central del límite, es un montaje experimental factible y regular. ¿Y qué condiciones de frontera?

Si consideramos de Dirichlet, a una tira de longitud <math> 2a </math> le colocaremos en los extremos un material que se mantenga siempre a la misma temperatura, pongamos que a una temperatura nula. El sistema resulta

<math>\quad
\begin{align}
\begin{cases} u_{t} - u_{xx} = 0 \quad &\forall t > 0, x\in [-a,a],\\
u(-a,t) = u(a,t) = 0 \quad &\forall t > 0, \\
u(x,0) = e^{-x^2} \quad &\forall x\in [-a,a].
\end{cases}
\end{align}
</math>

La solución al resolver el problema de autofunciones viene dada por

<math>\quad
u(x,t) = \sum_{n=1}^\infty A_{n} \cos\left(\frac{\pi}{a}(n - \frac{1}{2}) x\right) e^{-\dfrac{\pi^2}{a^2}\left(n - \frac{1}{2}\right)^2 t},
\quad \text{ donde } \quad
A_n = \frac{\int_{-a}^{a} u(x,0) \cos\left(\frac{\pi}{a}(n - \frac{1}{2}) x\right)dx}{\int_{-a}^{a} \cos^2\left(\frac{\pi}{a}(n - \frac{1}{2}) x\right)dx}.
</math>

Consideremos ahora condiciones de frontera Neumann, con un flujo nulo que físicamente viene de poner extremos aislantes. Para cierto <math> b > 0,</math> y no necesariamente igual que con <math> a </math>. Consideraremos la ecuación del calor

<math>\quad
\begin{align}
\begin{cases} u_{t} - u_{xx} = 0 \quad &\forall t > 0 , x\in [-b,b],\\
u_x(-b,t) = u_x(b,t) = 0 \quad &\forall t > 0, \\
u(x,0) = e^{-x^2} \quad &\forall x\in[-b,b].
\end{cases}
\end{align}
</math>

De nuevo podemos resolver el sistema, por ser homogéneo, mediante separación de variables, obteniendo

<math>\quad
u(x,t) = \sum_{m=0}^\infty B_{m} cos\left(\frac{m\pi}{b} x\right) e^{-\dfrac{m^2\pi^2}{b^2} t},\quad \text{ donde } \quad
B_m = \frac{\int_{-b}^{b} u(x,0) \cos\left(\frac{m\pi}{b} x\right)dx}{\int_{-b}^{b} \cos^2\left(\frac{m\pi}{b}x\right)dx}.
</math>

Un último detalle antes de proceder a los errores es la obtención simbólica de la solución no acotada.

<math>\quad
\begin{align}
\begin{cases} u_{t} - u_{xx} = 0 \quad &\forall t>0, x\in \mathbb{R},\\
u(x,0) = e^{-x^2} \quad &\forall x\in \mathbb{R}
\end{cases}
\end{align}
\quad
\Longrightarrow
\quad
u(x,t) = \frac{1}{\sqrt{4\pi t}}e^{\frac{-x^2}{4t}} * e^{-x^2} = \frac{e^{\frac{-x^2}{4t+1}}}{\sqrt{4t+1}}.
</math>

Queremos comprobar cuán bien aproximan las soluciones del caso acotado Dirichlet y Neumann a la obtenida por la solución fundamental. Para ello, escojamos un tiempo máximo <math>t_{max} </math> hasta el cuál tomaremos el error en norma <math> L^{\infty} </math>. Queremos estudiar:

* El efecto de aumentar el rango de definición con <math> a,b </math>.
* El efecto de aumentar el número de términos en las series de los sistemas acotados.

=== Análisis de los errores ===
[[Archivo:GrupoADMR Calor Errores1.jpg|550px|thumb|right| Representación gráfica de las soluciones del problema acotado (<math>a = b = 20</math>) y no acotado truncando las series de Fourier en su término <math>100</math>, para un tiempo máximo de <math>1</math>; y de sus errores. Nótese el orden de magnitud de los errores -pequeño-, los picos generados por los errores numéricos dados por la fórmula del trapecio y la suavidad en la frontera dada por las condiciones en esta.]]

Consideramos los parámetros, donde <math> n </math> es el término de la serie de Fourier donde truncamos:

*<math> a = b = 20, </math>
*<math> t_{max} = 1, </math>
*<math> n = 100. </math>

Con <math> t_{max} </math> fijo podemos observar que al aumentar <math>a, b</math> y <math>n</math> el error entre las soluciones acotadas y no acotada disminuye. Nótese que se estanca a partir de cierto <math>n</math>. De hecho, el error llega a aumentar a partir de ciertos <math> a,b </math> si truncamos demasiado pronto o demasiado tarde por la aproximación del trapecio para particiones del espacio demasiado bastas -no observado por pericia en código final-.

{|style="margin: 0 auto;"
| [[Archivo:GrupoADMR Calor Errores MaxD.jpeg|250px|thumb|upright|Con <math>t_{max}=1</math>, vamos obteniendo errores máximos en escala logarítmica con diferentes longitudes y truncamientos entre la solución del problema no acotado y el acotado Dirichlet. Los errores son bastante pequeños si se toman suficientes términos de la serie y una discretización espacial adecuada (en otro caso, aunque no esté representado en esta versión final, los errores divergen)]]
| [[Archivo:GrupoADMR Calor Errores MaxNe.jpeg|250px|thumb|upright|Análogo a la figura adyacente. Nótese que ambas son buenas aproximaciones con suficiente longitud y términos de la serie.]]
|}

Al ampliar el tiempo con la <math>n</math> y fijada <math>a=b=15</math>, al hacer más grande <math>t_{max}</math> los errores máximos no aumentan. En cambio, si tomamos el máximo error para diferentes tiempos vemos con suficiente precisión que no hay diferencia pero que para truncamientos precoces Neumann tiene ligeramente menor error.

{|style="margin: 0 auto;"
| [[Archivo:GrupoADMR Calor Errores DNo.jpeg|250px|thumb|upright|Con dominio fijo, el error en escala logarítmica entre solución no acotada y sistema Dirichlet decae incluso con el paso del tiempo. Esto encaja con que el estado estacionario de Dirichlet será la constante de temperatura nula, al igual que la convolución con la solución fundamental.]]
| [[Archivo:GrupoADMR Calor Errores NoNe.jpeg|250px|thumb|upright|Con condiciones análogas en frontera Neumann, el resultado es igualmente excelente para aproximar el caso no acotado. Aunque en el caso Neumann la solución estacionaria no sea nula, sino un reparto equitativo del calor inicial, cuando el rango es suficientemente amplio, la aproximación será adecuada.]]
| [[Archivo:GrupoADMR Calor Errores DNe.jpeg|250px|thumb|upright|Podemos ver además que ambos sistemas se asemejan cada vez más, por lo que ambas aproximaciones son igualmente útiles para el caso no acotado.]]
|}

=Coeficiente de difusión=

[[Archivo:Fundamental.gif|300px|thumb|right|Solución fundamental de la ecuación del calor unidimensional para distintos valores del coeficiente de difusión. La representación gráfica parte del tiempo inicial <math/> 0.01, <math/> esperando una delta de Dirac en el instante inicial. Nótese que en este tiempo, a mayor coeficiente de difusión más se ha difundido el calor en el sistema, lo que se nota especialmente en el valor máximo observado y su disminución cuanto más crece el coeficiente.]]

Interpretemos el valor de <math/> D <math/> aprovechando el montaje físico anterior. Numéricamente vemos que cuanto mayor es este coeficiente más rápido se difunde el calor, evidenciado en la caída de la solución en tiempos cada vez más cortos.

[[Archivo:Convoluciones_teorica.gif|500px|thumb|right|Soluciones de la ecuación del calor unidimensional para distintos valores del coeficiente de difusión.]]

En el caso unidimensional acotado con solución

<math>
u(x,t) = \sum_{n=1}^{\infty} A_n e^{-D (n\pi)^2 t} \sin(n\pi x).
</math>

Podemos ver que un mayor coeficiente <math/> D <math/> acelera el decaimiento de la exponencial por lo que tendrá un efecto similar justificando su nombre.

= Códigos =

<source lang="matlab">
%%% Código generado junto con Chat GPT para la representación de la solución fundamental de la ecuación del calor unidimensional para distintos valores del coeficiente de difusión

clc
clear all
close all

% Vector de valores de D (coeficiente de difusión)
DD = [1,2,3,4,5,6,7,8,9,10];

% Nombre del archivo GIF de salida
output_gif = 'Fundamental.gif';

% Crear la figura
figura = figure(1);
grid on;
view(3);

% Bucle para iterar sobre cada valor de D
for i = 1:length(DD)
D = DD(i);

% Función fundamental de la ecuación del calor
phi_D = @(x, t) exp(-x.^2./(4*D*t))./sqrt(4*D*pi*t);

% Definición de la malla
xx = -1:0.01:1; % Rango de x
tt = 10^(-2):0.01:1; % Rango de t
[X, T] = meshgrid(xx, tt); % Malla de coordenadas
PHI_D = phi_D(X, T); % Evaluación de la función

% Graficar la superficie
clf;
surf(X, T, PHI_D);
shading flat
shading interp
colormap('jet')
xlabel('x');
ylabel('t');
zlabel('$\Phi_D(x,t)$', 'Interpreter', 'latex');
title(['Solución fundamental de la ecuación del calor, cuando D=', num2str(D)])
axis([-1, 1, 0, 1, 0, 3])
drawnow;

% Captura de la imagen para el GIF
frame = getframe(figura);
img = frame2im(frame);
[imind, cm] = rgb2ind(img, 256);

% Escribir la imagen en el archivo GIF
if i == 1
imwrite(imind, cm, output_gif, 'gif', 'LoopCount', Inf, 'DelayTime', 0.7);
else
imwrite(imind, cm, output_gif, 'gif', 'WriteMode', 'append', 'DelayTime', 0.7);
end

end
</source>

<source lang="matlab">
%%% Código generado junto con Chat GPT para la representación de la soluciones de la ecuación del calor unidimensional para distintos valores del coeficiente de difusión

clc
clear all
close all

% Vector de valores de D (coeficiente de difusión)
DD = [1,10,20,30,40,50,60,70,80,90,100];

% Función gaussiana
gaus = @(x) (exp(-x.^2));

% Nombre del archivo GIF de salida
output_gif = 'Convoluciones_teorica.gif';

% Crear la figura
figura = figure(1);
grid on;
view(3);

% Bucle para iterar sobre cada valor de D
for h = 1:length(DD)
D = DD(h);

% Funciones fundamentales y convolución
phi_D = @(x, t) exp(-x.^2./(4*D*t))./sqrt(4*D*pi*t);
convo = @(x,t) (exp(-x.^2./(4*D.*t+1))/sqrt(4*D.*t+1));

% Definición de la malla
xx = -10:0.1:10; % Rango de x
tt = 10^(-2):0.01:1; % Rango de t

% Inicialización de la matriz de convolución
conv = zeros(length(xx), length(tt));

% Cálculo de la convolución
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
val_conv(i,j) = convo(auxx, auxt);
end
end

% Graficar la superficie
clf;
hold on;
view(3);
surf(xx, tt, val_conv.');
xlabel('x');
ylabel('t');
zlabel('$\Phi_D(x,t)*N(x)$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap('jet');
title(['Soluciones en dominio no acotado, cuando D=', num2str(D)]);
axis([-10, 10, 0, 1, 0, 1]);
hold off;
drawnow;

% Captura de la imagen para el GIF
frame = getframe(figura);
img = frame2im(frame);
[imind, cm] = rgb2ind(img, 256);

% Escribir la imagen en el archivo GIF
if h == 1
imwrite(imind, cm, output_gif, 'gif', 'LoopCount', Inf, 'DelayTime', 0.7);
else
imwrite(imind, cm, output_gif, 'gif', 'WriteMode', 'append', 'DelayTime', 0.7);
end

end
</source>

<source lang="matlab">
%%% Código generado junto con Chat GPT para la representación de los errores entre las soluciones

function [Diff12,Difft1,Difft2] = errores2(n_max, a, b,t_max)

% Función como condición inicial: una función Gaussiana
gaus = @(x) (exp(-x.^2));
u0 = @(x) exp(-x.^2);
D = 1; % Coeficiente de difusión

% Funciones fundamentales y convolución para el modelo de difusión
phi_D = @(x, t) exp(-x.^2./(4*D*t))./sqrt(4*D*pi*t);
convo = @(x,t) (exp(-x.^2./(4*D.*t+1))/sqrt(4*D.*t+1));

% Definición de la malla (espacio y tiempo)
xx = -a:0.01:a;
tt = 10^(-2):0.01:t_max;
[X, T] = meshgrid(xx, tt);

% Inicialización de la matriz de convolución
conv = zeros(length(xx), length(tt));

% Cálculo de la convolución
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
val_conv(j,i) = convo(auxx, auxt);
end
end

% Inicialización de la matriz de resultados para la solución Dirichlet
res_dir = zeros(length(xx), length(tt));

% Cálculo de la solución directa
for n=1:n_max
A(n) = trapz(xx, u0(xx) .* cos((pi/a)*(n - 1/2) * xx)) /...
trapz(xx, cos((pi/a)*(n - 1/2) * xx).^2);
end

% Evaluación de la solución directa en el dominio espacio-temporal
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
u = 0;

% Sumar los términos de la serie de Fourier
for n=1:n_max
u = u + A(n) .* cos((pi/a)*(n - 1/2) * auxx) .*...
exp(-((pi^2 / a^2) * (n - 1/2)^2) .* auxt);
end
res_dir(i,j) = u;
end
end

% Inicialización de la matriz de resultados para la solución en frontera Neumann
res_neu = zeros(length(xx), length(tt));

% Cálculo de la solución en frontera Neumann
for n=1:n_max
B(n) = trapz(xx, u0(xx) .* cos((n*pi/b) * xx)) /...
trapz(xx, cos((n*pi/b) * xx).^2);
end

% Evaluación de la solución en frontera Neumann en el dominio espacio-temporal
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
u = trapz(xx, u0(xx))./(2.*b);

% Sumar los términos de la serie de Fourier
for n=1:n_max
u = u + B(n) .* cos((n*pi/b) * auxx) .* exp(-((n^2 * pi^2 / b^2)) .* auxt);
end
res_neu(i,j) = u;
end
end

% Transponer los resultados para compararlos con la convolución
U1 = res_dir.';
U2 = res_neu.';

% Calcular las diferencias entre las soluciones y la convolución
Diff12 = abs(U1 - U2);
Difft1 = abs(val_conv - U1);
Difft2 = abs(val_conv - U2);

end
</source>

<source lang="matlab">
%%% Código generado junto con Chat GPT para la representación de las soluciones de la ecuación del calor unidimensional junto con la comparación de sus errores

clc
clear all
close all

% Función como condicion inicial
% Se define la función inicial gausiana y la función u0 que se utilizarán más adelante para los cálculos.
gaus = @(x) (exp(-x.^2));
u0 = @(x) exp(-x.^2);
D = 1;

% Funciones fundamentales y convolución
% Se definen las funciones que describen la solución fundamental de la ecuación y la convolución utilizada.
phi_D = @(x, t) exp(-x.^2./(4*D*t))./sqrt(4*D*pi*t);
convo = @(x,t) (exp(-x.^2./(4*D.*t+1))/sqrt(4*D.*t+1));

% Definición de la malla
% Se definen los parámetros de la malla para las variables espaciales (xx) y temporales (tt), y se crea la malla 2D.
a = 100;
b = 100;
tmax = 1;
n_max = 1000;

xx = linspace(-a, a, 1000);
tt = linspace(10^(-2), tmax, 1000);
[X, T] = meshgrid(xx, tt);

% Inicialización de la matriz de convolución
% Se crea una matriz vacía para almacenar los valores de la convolución.
conv = zeros(length(xx), length(tt));

% Cálculo de la convolución
% Se calcula la convolución para cada punto en la malla de la variable temporal tt y espacial xx.
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
val_conv(i,j) = convo(auxx, auxt);
end
end

% Resolución de la ecuación con condiciones de frontera Dirichlet
% Se calculan los coeficientes A para la expansión en series de Fourier y se calcula la solución directa en cada paso temporal.
res_dir = zeros(length(xx), length(tt));

for n=1:n_max
A(n) = trapz(xx, u0(xx) .* cos((pi/a)*(n - 1/2) * xx)) /...
trapz(xx, cos((pi/a)*(n - 1/2) * xx).^2);
end

for j = 1:length(tt)
auxt = tt(j);

for i = 1:length(xx)
auxx = xx(i);
u = 0;

for n=1:n_max
u = u + A(n) .* cos((pi/a)*(n - 1/2) * auxx) .*...
exp(-((pi^2 / a^2) * (n - 1/2)^2) .* auxt);

end
res_dir(j,i) = u;
end
end

% Resolución de la ecuación con condiciones de frontera Neumann
% Se calculan los coeficientes B para la expansión en series de Fourier con las condiciones de frontera Neumann y se calcula la solución en cada paso temporal.
res_neu = zeros(length(xx), length(tt));

for n=1:n_max
B(n) = trapz(xx, u0(xx) .* cos((n*pi/b) * xx)) /...
trapz(xx, cos((n*pi/b) * xx).^2);
end

for j = 1:length(tt)
auxt = tt(j);

for i = 1:length(xx)
auxx = xx(i);
u = trapz(xx, u0(xx))./(2.*b);

for n=1:n_max
u = u + B(n) .* cos((n*pi/b) * auxx) .* exp(-((n^2 * pi^2 / b^2)) .* auxt);

end
res_neu(j,i) = u;
end
end

% Comparación entre las soluciones obtenidas
% Se calcula la diferencia entre las soluciones directa y Neumann, así como entre las soluciones y la convolución.
U1 = res_dir;
U2 = res_neu;

Diff12 = abs(U1 - U2);
Difft1 = abs(val_conv.' - U1);
Difft2 = abs(val_conv.' - U2);

% Visualización de los resultados
% Se muestran las superficies 3D de las soluciones y sus diferencias en gráficos para comparar los resultados obtenidos.
figure(1);

subplot(2,4,2);
view(3);
surf(X, T, U1);
xlabel('x');
ylabel('t');
zlabel('$u_1(x,t)$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca,'jet');
axis([-a,a,0,tmax,0,1])
title('$u_1(x,t)$', 'Interpreter', 'latex');

subplot(2,4,3);
view(3);
surf(X, T, U2);
xlabel('x');
ylabel('t');
zlabel('$u_2(x,t)$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca, 'jet');
axis([-b,b,0,tmax,0,1])
title('$u_2(x,t)$', 'Interpreter', 'latex');

subplot(2,4,4);
view(3);
surf(X, T, Diff12);
xlabel('x');
ylabel('t');
zlabel('$|u_1(x,t)-u_2(x,t)|$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca, 'hot');
title('$|u_1(x,t)-u_2(x,t)|$', 'Interpreter', 'latex');

subplot(2,4,5)
view(3);
surf(X, T, val_conv.');
xlabel('x');
ylabel('t');
zlabel('$\Phi_D(x,t)*N(x)$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca, 'jet');
axis([-a,a,0,tmax,0,1])
title('$\Phi_D(x,t)*N(x)$', 'Interpreter', 'latex');

subplot(2,4,6);
view(3);
surf(X, T, Difft1);
xlabel('x');
ylabel('t');
zlabel('$|\Phi_D(x,t)*N(x)-u_1(x,t)|$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca, 'hot');
title('$|\Phi_D(x,t)*N(x)-u_1(x,t)|$', 'Interpreter', 'latex');

subplot(2,4,7);
view(3);
surf(X, T, Difft2);
xlabel('x');
ylabel('t');
zlabel('$|\Phi_D(x,t)*N(x)-u_2(x,t)|$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca, 'hot');
title('$|\Phi_D(x,t)*N(x)-u_2(x,t)|$', 'Interpreter', 'latex');
</source>

= Referencias =

* Sandro Salsa, Gianmaria Verzini. Partial Differential Equations in Action From Modelling to Theory ([https://link.springer.com/book/10.1007/978-3-031-21853-8 [1<nowiki>]</nowiki>]).

[[Categoría:EDP]]
[[Categoría:EDP24/25]]

Ecuación del calor (ADMR)

2025-03-19T07:38:35Z

Daniel Rodríguez Calderón: /* Códigos */

{{ TrabajoED | Ecuación del calor (Grupo ADMR). | [[:Categoría:EDP|EDP]]|[[:Categoría:EDP24/25|2024-25]] |
*Ángel De Lucas Miranda,
*Daniel Rodríguez Calderón,
*Marcos Cabellos Hernández,
*Rafael Pascual Ortega.}}

=Introducción=

El desarrollo de la ecuación del calor se atribuye principalmente al matemático y físico francés Joseph Fourier, quien en 1822 publicó su obra ''Théorie analytique de la chaleur''. En este trabajo formuló una ecuación diferencial en derivadas parciales que describe la propagación del calor por difusión en un medio continuo. En esta ocasión, profundicemos en dos aspectos de este modelo.

Desde un punto de vista físico, estas ecuaciones permiten la modelización de fenómenos como el calentamiento de materiales o la difusión de solutos en fluidos estacionarios. Sin embargo, la influencia de la física es más profunda: de hecho, el coeficiente de difusión lleva su nombre por la interpretación física de su acción en el sistema. ¿Cuál es este efecto, y por qué se lo denota de esta forma?

Planteemos también el paso entre problemas de dominio acotado y los de no acotado. ¿Podemos aproximar las soluciones del sistema no acotado mediante las del acotado? Y de ser así, ¿qué condiciones de frontera en una dimensión funcionan mejor?

=Solución acotada vs Solución no acotada=
Demos por conocidos el uso de separación de variables para resolver problemas de autofunciones en caso acotado y el manejo de la solución fundamental junto con la convolución para el sistema no acotado. Para comparar todas las soluciones correspondientes, consideremos sistema <math> 1-</math>dimensional, con una gaussiana <math> e^{-x^2} </math> como condición inicial. Justificamos físicamente como el uso de un láser para calentar un punto de una tira fina de metal. Considerando que los fotones puedan desviarse como una distribución normal por el teorema central del límite, es un montaje experimental factible y regular. ¿Y qué condiciones de frontera?

Si consideramos de Dirichlet, a una tira de longitud <math> 2a </math> le colocaremos en los extremos un material que se mantenga siempre a la misma temperatura, pongamos que a una temperatura nula. El sistema resulta

<math>\quad
\begin{align}
\begin{cases} u_{t} - u_{xx} = 0 \quad &\forall t > 0, x\in [-a,a],\\
u(-a,t) = u(a,t) = 0 \quad &\forall t > 0, \\
u(x,0) = e^{-x^2} \quad &\forall x\in [-a,a].
\end{cases}
\end{align}
</math>

La solución al resolver el problema de autofunciones viene dada por

<math>\quad
u(x,t) = \sum_{n=1}^\infty A_{n} \cos\left(\frac{\pi}{a}(n - \frac{1}{2}) x\right) e^{-\dfrac{\pi^2}{a^2}\left(n - \frac{1}{2}\right)^2 t},
\quad \text{ donde } \quad
A_n = \frac{\int_{-a}^{a} u(x,0) \cos\left(\frac{\pi}{a}(n - \frac{1}{2}) x\right)dx}{\int_{-a}^{a} \cos^2\left(\frac{\pi}{a}(n - \frac{1}{2}) x\right)dx}.
</math>

Consideremos ahora condiciones de frontera Neumann, con un flujo nulo que físicamente viene de poner extremos aislantes. Para cierto <math> b > 0,</math> y no necesariamente igual que con <math> a </math>. Consideraremos la ecuación del calor

<math>\quad
\begin{align}
\begin{cases} u_{t} - u_{xx} = 0 \quad &\forall t > 0 , x\in [-b,b],\\
u_x(-b,t) = u_x(b,t) = 0 \quad &\forall t > 0, \\
u(x,0) = e^{-x^2} \quad &\forall x\in[-b,b].
\end{cases}
\end{align}
</math>

De nuevo podemos resolver el sistema, por ser homogéneo, mediante separación de variables, obteniendo

<math>\quad
u(x,t) = \sum_{m=0}^\infty B_{m} cos\left(\frac{m\pi}{b} x\right) e^{-\dfrac{m^2\pi^2}{b^2} t},\quad \text{ donde } \quad
B_m = \frac{\int_{-b}^{b} u(x,0) \cos\left(\frac{m\pi}{b} x\right)dx}{\int_{-b}^{b} \cos^2\left(\frac{m\pi}{b}x\right)dx}.
</math>

Un último detalle antes de proceder a los errores es la obtención simbólica de la solución no acotada.

<math>\quad
\begin{align}
\begin{cases} u_{t} - u_{xx} = 0 \quad &\forall t>0, x\in \mathbb{R},\\
u(x,0) = e^{-x^2} \quad &\forall x\in \mathbb{R}
\end{cases}
\end{align}
\quad
\Longrightarrow
\quad
u(x,t) = \frac{1}{\sqrt{4\pi t}}e^{\frac{-x^2}{4t}} * e^{-x^2} = \frac{e^{\frac{-x^2}{4t+1}}}{\sqrt{4t+1}}.
</math>

Queremos comprobar cuán bien aproximan las soluciones del caso acotado Dirichlet y Neumann a la obtenida por la solución fundamental. Para ello, escojamos un tiempo máximo <math>t_{max} </math> hasta el cuál tomaremos el error en norma <math> L^{\infty} </math>. Queremos estudiar:

* El efecto de aumentar el rango de definición con <math> a,b </math>.
* El efecto de aumentar el número de términos en las series de los sistemas acotados.

=== Análisis de los errores ===
[[Archivo:GrupoADMR Calor Errores1.jpg|550px|thumb|right| Representación gráfica de las soluciones del problema acotado (<math>a = b = 20</math>) y no acotado truncando las series de Fourier en su término <math>100</math>, para un tiempo máximo de <math>1</math>; y de sus errores. Nótese el orden de magnitud de los errores -pequeño-, los picos generados por los errores numéricos dados por la fórmula del trapecio y la suavidad en la frontera dada por las condiciones en esta.]]

Consideramos los parámetros, donde <math> n </math> es el término de la serie de Fourier donde truncamos:

*<math> a = b = 20, </math>
*<math> t_{max} = 1, </math>
*<math> n = 100. </math>

Con <math> t_{max} </math> fijo podemos observar que al aumentar <math>a, b</math> y <math>n</math> el error entre las soluciones acotadas y no acotada disminuye. Nótese que se estanca a partir de cierto <math>n</math>. De hecho, el error llega a aumentar a partir de ciertos <math> a,b </math> si truncamos demasiado pronto o demasiado tarde por la aproximación del trapecio para particiones del espacio demasiado bastas -no observado por pericia en código final-.

{|style="margin: 0 auto;"
| [[Archivo:GrupoADMR Calor Errores MaxD.jpeg|250px|thumb|upright|Con <math>t_{max}=1</math>, vamos obteniendo errores máximos en escala logarítmica con diferentes longitudes y truncamientos entre la solución del problema no acotado y el acotado Dirichlet. Los errores son bastante pequeños si se toman suficientes términos de la serie y una discretización espacial adecuada (en otro caso, aunque no esté representado en esta versión final, los errores divergen)]]
| [[Archivo:GrupoADMR Calor Errores MaxNe.jpeg|250px|thumb|upright|Análogo a la figura adyacente. Nótese que ambas son buenas aproximaciones con suficiente longitud y términos de la serie.]]
|}

Al ampliar el tiempo con la <math>n</math> y fijada <math>a=b=15</math>, al hacer más grande <math>t_{max}</math> los errores máximos no aumentan. En cambio, si tomamos el máximo error para diferentes tiempos vemos con suficiente precisión que no hay diferencia pero que para truncamientos precoces Neumann tiene ligeramente menor error.

{|style="margin: 0 auto;"
| [[Archivo:GrupoADMR Calor Errores DNo.jpeg|250px|thumb|upright|Con dominio fijo, el error en escala logarítmica entre solución no acotada y sistema Dirichlet decae incluso con el paso del tiempo. Esto encaja con que el estado estacionario de Dirichlet será la constante de temperatura nula, al igual que la convolución con la solución fundamental.]]
| [[Archivo:GrupoADMR Calor Errores NoNe.jpeg|250px|thumb|upright|Con condiciones análogas en frontera Neumann, el resultado es igualmente excelente para aproximar el caso no acotado. Aunque en el caso Neumann la solución estacionaria no sea nula, sino un reparto equitativo del calor inicial, cuando el rango es suficientemente amplio, la aproximación será adecuada.]]
| [[Archivo:GrupoADMR Calor Errores DNe.jpeg|250px|thumb|upright|Podemos ver además que ambos sistemas se asemejan cada vez más, por lo que ambas aproximaciones son igualmente útiles para el caso no acotado.]]
|}

=Coeficiente de difusión=

[[Archivo:Fundamental.gif|300px|thumb|right|Solución fundamental de la ecuación del calor unidimensional para distintos valores del coeficiente de difusión. La representación gráfica parte del tiempo inicial <math/> 0.01, <math/> esperando una delta de Dirac en el instante inicial. Nótese que en este tiempo, a mayor coeficiente de difusión más se ha difundido el calor en el sistema, lo que se nota especialmente en el valor máximo observado y su disminución cuanto más crece el coeficiente.]]

Interpretemos el valor de <math/> D <math/> aprovechando el montaje físico anterior. Numéricamente vemos que cuanto mayor es este coeficiente más rápido se difunde el calor, evidenciado en la caída de la solución en tiempos cada vez más cortos.

[[Archivo:Convoluciones_teorica.gif|500px|thumb|right|Soluciones de la ecuación del calor unidimensional para distintos valores del coeficiente de difusión.]]

En el caso unidimensional acotado con solución

<math>
u(x,t) = \sum_{n=1}^{\infty} A_n e^{-D (n\pi)^2 t} \sin(n\pi x).
</math>

Podemos ver que un mayor coeficiente <math/> D <math/> acelera el decaimiento de la exponencial por lo que tendrá un efecto similar justificando su nombre.

= Códigos =

<source lang="matlab">
%%% Código generado junto con Chat GPT para la representación de la solución fundamental de la ecuación del calor unidimensional para distintos valores del coeficiente de difusión

clc
clear all
close all

% Vector de valores de D (coeficiente de difusión)
DD = [1,2,3,4,5,6,7,8,9,10];

% Nombre del archivo GIF de salida
output_gif = 'Fundamental.gif';

% Crear la figura
figura = figure(1);
grid on;
view(3);

% Bucle para iterar sobre cada valor de D
for i = 1:length(DD)
D = DD(i);

% Función fundamental de la ecuación del calor
phi_D = @(x, t) exp(-x.^2./(4*D*t))./sqrt(4*D*pi*t);

% Definición de la malla
xx = -1:0.01:1; % Rango de x
tt = 10^(-2):0.01:1; % Rango de t
[X, T] = meshgrid(xx, tt); % Malla de coordenadas
PHI_D = phi_D(X, T); % Evaluación de la función

% Graficar la superficie
clf;
surf(X, T, PHI_D);
shading flat
shading interp
colormap('jet')
xlabel('x');
ylabel('t');
zlabel('$\Phi_D(x,t)$', 'Interpreter', 'latex');
title(['Solución fundamental de la ecuación del calor, cuando D=', num2str(D)])
axis([-1, 1, 0, 1, 0, 3])
drawnow;

% Captura de la imagen para el GIF
frame = getframe(figura);
img = frame2im(frame);
[imind, cm] = rgb2ind(img, 256);

% Escribir la imagen en el archivo GIF
if i == 1
imwrite(imind, cm, output_gif, 'gif', 'LoopCount', Inf, 'DelayTime', 0.7);
else
imwrite(imind, cm, output_gif, 'gif', 'WriteMode', 'append', 'DelayTime', 0.7);
end

end
</source>

<source lang="matlab">
%%% Código generado junto con Chat GPT para la representación de la soluciones de la ecuación del calor unidimensional para distintos valores del coeficiente de difusión

clc
clear all
close all

% Vector de valores de D (coeficiente de difusión)
DD = [1,10,20,30,40,50,60,70,80,90,100];

% Función gaussiana
gaus = @(x) (exp(-x.^2));

% Nombre del archivo GIF de salida
output_gif = 'Convoluciones_teorica.gif';

% Crear la figura
figura = figure(1);
grid on;
view(3);

% Bucle para iterar sobre cada valor de D
for h = 1:length(DD)
D = DD(h);

% Funciones fundamentales y convolución
phi_D = @(x, t) exp(-x.^2./(4*D*t))./sqrt(4*D*pi*t);
convo = @(x,t) (exp(-x.^2./(4*D.*t+1))/sqrt(4*D.*t+1));

% Definición de la malla
xx = -10:0.1:10; % Rango de x
tt = 10^(-2):0.01:1; % Rango de t

% Inicialización de la matriz de convolución
conv = zeros(length(xx), length(tt));

% Cálculo de la convolución
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
val_conv(i,j) = convo(auxx, auxt);
end
end

% Graficar la superficie
clf;
hold on;
view(3);
surf(xx, tt, val_conv.');
xlabel('x');
ylabel('t');
zlabel('$\Phi_D(x,t)*N(x)$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap('jet');
title(['Soluciones en dominio no acotado, cuando D=', num2str(D)]);
axis([-10, 10, 0, 1, 0, 1]);
hold off;
drawnow;

% Captura de la imagen para el GIF
frame = getframe(figura);
img = frame2im(frame);
[imind, cm] = rgb2ind(img, 256);

% Escribir la imagen en el archivo GIF
if h == 1
imwrite(imind, cm, output_gif, 'gif', 'LoopCount', Inf, 'DelayTime', 0.7);
else
imwrite(imind, cm, output_gif, 'gif', 'WriteMode', 'append', 'DelayTime', 0.7);
end

end
</source>

<source lang="matlab">
%%% Código generado junto con Chat GPT para la representación de los errores entre las soluciones

function [Diff12,Difft1,Difft2] = errores2(n_max, a, b,t_max)

% Función como condición inicial: una función Gaussiana
gaus = @(x) (exp(-x.^2));
u0 = @(x) exp(-x.^2);
D = 1; % Coeficiente de difusión

% Funciones fundamentales y convolución para el modelo de difusión
phi_D = @(x, t) exp(-x.^2./(4*D*t))./sqrt(4*D*pi*t);
convo = @(x,t) (exp(-x.^2./(4*D.*t+1))/sqrt(4*D.*t+1));

% Definición de la malla (espacio y tiempo)
xx = -a:0.01:a;
tt = 10^(-2):0.01:t_max;
[X, T] = meshgrid(xx, tt);

% Inicialización de la matriz de convolución
conv = zeros(length(xx), length(tt));

% Cálculo de la convolución
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
val_conv(j,i) = convo(auxx, auxt);
end
end

% Inicialización de la matriz de resultados para la solución Dirichlet
res_dir = zeros(length(xx), length(tt));

% Cálculo de la solución directa
for n=1:n_max
A(n) = trapz(xx, u0(xx) .* cos((pi/a)*(n - 1/2) * xx)) /...
trapz(xx, cos((pi/a)*(n - 1/2) * xx).^2);
end

% Evaluación de la solución directa en el dominio espacio-temporal
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
u = 0;

% Sumar los términos de la serie de Fourier
for n=1:n_max
u = u + A(n) .* cos((pi/a)*(n - 1/2) * auxx) .*...
exp(-((pi^2 / a^2) * (n - 1/2)^2) .* auxt);
end
res_dir(i,j) = u;
end
end

% Inicialización de la matriz de resultados para la solución en frontera Neumann
res_neu = zeros(length(xx), length(tt));

% Cálculo de la solución en frontera Neumann
for n=1:n_max
B(n) = trapz(xx, u0(xx) .* cos((n*pi/b) * xx)) /...
trapz(xx, cos((n*pi/b) * xx).^2);
end

% Evaluación de la solución en frontera Neumann en el dominio espacio-temporal
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
u = trapz(xx, u0(xx))./(2.*b);

% Sumar los términos de la serie de Fourier
for n=1:n_max
u = u + B(n) .* cos((n*pi/b) * auxx) .* exp(-((n^2 * pi^2 / b^2)) .* auxt);
end
res_neu(i,j) = u;
end
end

% Transponer los resultados para compararlos con la convolución
U1 = res_dir.';
U2 = res_neu.';

% Calcular las diferencias entre las soluciones y la convolución
Diff12 = abs(U1 - U2);
Difft1 = abs(val_conv - U1);
Difft2 = abs(val_conv - U2);

end
</source>

<source lang="matlab">
%%% Código generado junto con Chat GPT para la representación de las soluciones de la ecuación del calor unidimensional junto con la comparación de sus errores

clc
clear all
close all

% Función como condicion inicial
% Se define la función inicial gausiana y la función u0 que se utilizarán más adelante para los cálculos.
gaus = @(x) (exp(-x.^2));
u0 = @(x) exp(-x.^2);
D = 1;

% Funciones fundamentales y convolución
% Se definen las funciones que describen la solución fundamental de la ecuación y la convolución utilizada.
phi_D = @(x, t) exp(-x.^2./(4*D*t))./sqrt(4*D*pi*t);
convo = @(x,t) (exp(-x.^2./(4*D.*t+1))/sqrt(4*D.*t+1));

% Definición de la malla
% Se definen los parámetros de la malla para las variables espaciales (xx) y temporales (tt), y se crea la malla 2D.
a = 100;
b = 100;
tmax = 1;
n_max = 1000;

xx = linspace(-a, a, 1000);
tt = linspace(10^(-2), tmax, 1000);
[X, T] = meshgrid(xx, tt);

% Inicialización de la matriz de convolución
% Se crea una matriz vacía para almacenar los valores de la convolución.
conv = zeros(length(xx), length(tt));

% Cálculo de la convolución
% Se calcula la convolución para cada punto en la malla de la variable temporal tt y espacial xx.
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
val_conv(i,j) = convo(auxx, auxt);
end
end

% Resolución de la ecuación con condiciones de frontera Dirichlet
% Se calculan los coeficientes A para la expansión en series de Fourier y se calcula la solución directa en cada paso temporal.
res_dir = zeros(length(xx), length(tt));

for n=1:n_max
A(n) = trapz(xx, u0(xx) .* cos((pi/a)*(n - 1/2) * xx)) /...
trapz(xx, cos((pi/a)*(n - 1/2) * xx).^2);
end

for j = 1:length(tt)
auxt = tt(j);

for i = 1:length(xx)
auxx = xx(i);
u = 0;

for n=1:n_max
u = u + A(n) .* cos((pi/a)*(n - 1/2) * auxx) .*...
exp(-((pi^2 / a^2) * (n - 1/2)^2) .* auxt);

end
res_dir(j,i) = u;
end
end

% Resolución de la ecuación con condiciones de frontera Neumann
% Se calculan los coeficientes B para la expansión en series de Fourier con las condiciones de frontera Neumann y se calcula la solución en cada paso temporal.
res_neu = zeros(length(xx), length(tt));

for n=1:n_max
B(n) = trapz(xx, u0(xx) .* cos((n*pi/b) * xx)) /...
trapz(xx, cos((n*pi/b) * xx).^2);
end

for j = 1:length(tt)
auxt = tt(j);

for i = 1:length(xx)
auxx = xx(i);
u = trapz(xx, u0(xx))./(2.*b);

for n=1:n_max
u = u + B(n) .* cos((n*pi/b) * auxx) .* exp(-((n^2 * pi^2 / b^2)) .* auxt);

end
res_neu(j,i) = u;
end
end

% Comparación entre las soluciones obtenidas
% Se calcula la diferencia entre las soluciones directa y Neumann, así como entre las soluciones y la convolución.
U1 = res_dir;
U2 = res_neu;

Diff12 = abs(U1 - U2);
Difft1 = abs(val_conv.' - U1);
Difft2 = abs(val_conv.' - U2);

% Visualización de los resultados
% Se muestran las superficies 3D de las soluciones y sus diferencias en gráficos para comparar los resultados obtenidos.
figure(1);

subplot(2,4,2);
view(3);
surf(X, T, U1);
xlabel('x');
ylabel('t');
zlabel('$u_1(x,t)$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca,'jet');
axis([-a,a,0,tmax,0,1])
title('$u_1(x,t)$', 'Interpreter', 'latex');

subplot(2,4,3);
view(3);
surf(X, T, U2);
xlabel('x');
ylabel('t');
zlabel('$u_2(x,t)$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca, 'jet');
axis([-b,b,0,tmax,0,1])
title('$u_2(x,t)$', 'Interpreter', 'latex');

subplot(2,4,4);
view(3);
surf(X, T, Diff12);
xlabel('x');
ylabel('t');
zlabel('$|u_1(x,t)-u_2(x,t)|$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca, 'hot');
title('$|u_1(x,t)-u_2(x,t)|$', 'Interpreter', 'latex');

subplot(2,4,5)
view(3);
surf(X, T, val_conv.');
xlabel('x');
ylabel('t');
zlabel('$\Phi_D(x,t)*N(x)$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca, 'jet');
axis([-a,a,0,tmax,0,1])
title('$\Phi_D(x,t)*N(x)$', 'Interpreter', 'latex');

subplot(2,4,6);
view(3);
surf(X, T, Difft1);
xlabel('x');
ylabel('t');
zlabel('$|\Phi_D(x,t)*N(x)-u_1(x,t)|$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca, 'hot');
title('$|\Phi_D(x,t)*N(x)-u_1(x,t)|$', 'Interpreter', 'latex');

subplot(2,4,7);
view(3);
surf(X, T, Difft2);
xlabel('x');
ylabel('t');
zlabel('$|\Phi_D(x,t)*N(x)-u_2(x,t)|$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca, 'hot');
title('$|\Phi_D(x,t)*N(x)-u_2(x,t)|$', 'Interpreter', 'latex');
</source>

= Referencias =

* Sandro Salsa, Gianmaria Verzini. Partial Differential Equations in Action From Modelling to Theory ([https://link.springer.com/book/10.1007/978-3-031-21853-8 [1<nowiki>]</nowiki>]).

[[Categoría:EDP]]
[[Categoría:EDP24/25]]

Ecuación del calor (ADMR)

2025-03-19T07:37:36Z

Daniel Rodríguez Calderón: /* Coeficiente de difusión */

{{ TrabajoED | Ecuación del calor (Grupo ADMR). | [[:Categoría:EDP|EDP]]|[[:Categoría:EDP24/25|2024-25]] |
*Ángel De Lucas Miranda,
*Daniel Rodríguez Calderón,
*Marcos Cabellos Hernández,
*Rafael Pascual Ortega.}}

=Introducción=

El desarrollo de la ecuación del calor se atribuye principalmente al matemático y físico francés Joseph Fourier, quien en 1822 publicó su obra ''Théorie analytique de la chaleur''. En este trabajo formuló una ecuación diferencial en derivadas parciales que describe la propagación del calor por difusión en un medio continuo. En esta ocasión, profundicemos en dos aspectos de este modelo.

Desde un punto de vista físico, estas ecuaciones permiten la modelización de fenómenos como el calentamiento de materiales o la difusión de solutos en fluidos estacionarios. Sin embargo, la influencia de la física es más profunda: de hecho, el coeficiente de difusión lleva su nombre por la interpretación física de su acción en el sistema. ¿Cuál es este efecto, y por qué se lo denota de esta forma?

Planteemos también el paso entre problemas de dominio acotado y los de no acotado. ¿Podemos aproximar las soluciones del sistema no acotado mediante las del acotado? Y de ser así, ¿qué condiciones de frontera en una dimensión funcionan mejor?

=Solución acotada vs Solución no acotada=
Demos por conocidos el uso de separación de variables para resolver problemas de autofunciones en caso acotado y el manejo de la solución fundamental junto con la convolución para el sistema no acotado. Para comparar todas las soluciones correspondientes, consideremos sistema <math> 1-</math>dimensional, con una gaussiana <math> e^{-x^2} </math> como condición inicial. Justificamos físicamente como el uso de un láser para calentar un punto de una tira fina de metal. Considerando que los fotones puedan desviarse como una distribución normal por el teorema central del límite, es un montaje experimental factible y regular. ¿Y qué condiciones de frontera?

Si consideramos de Dirichlet, a una tira de longitud <math> 2a </math> le colocaremos en los extremos un material que se mantenga siempre a la misma temperatura, pongamos que a una temperatura nula. El sistema resulta

<math>\quad
\begin{align}
\begin{cases} u_{t} - u_{xx} = 0 \quad &\forall t > 0, x\in [-a,a],\\
u(-a,t) = u(a,t) = 0 \quad &\forall t > 0, \\
u(x,0) = e^{-x^2} \quad &\forall x\in [-a,a].
\end{cases}
\end{align}
</math>

La solución al resolver el problema de autofunciones viene dada por

<math>\quad
u(x,t) = \sum_{n=1}^\infty A_{n} \cos\left(\frac{\pi}{a}(n - \frac{1}{2}) x\right) e^{-\dfrac{\pi^2}{a^2}\left(n - \frac{1}{2}\right)^2 t},
\quad \text{ donde } \quad
A_n = \frac{\int_{-a}^{a} u(x,0) \cos\left(\frac{\pi}{a}(n - \frac{1}{2}) x\right)dx}{\int_{-a}^{a} \cos^2\left(\frac{\pi}{a}(n - \frac{1}{2}) x\right)dx}.
</math>

Consideremos ahora condiciones de frontera Neumann, con un flujo nulo que físicamente viene de poner extremos aislantes. Para cierto <math> b > 0,</math> y no necesariamente igual que con <math> a </math>. Consideraremos la ecuación del calor

<math>\quad
\begin{align}
\begin{cases} u_{t} - u_{xx} = 0 \quad &\forall t > 0 , x\in [-b,b],\\
u_x(-b,t) = u_x(b,t) = 0 \quad &\forall t > 0, \\
u(x,0) = e^{-x^2} \quad &\forall x\in[-b,b].
\end{cases}
\end{align}
</math>

De nuevo podemos resolver el sistema, por ser homogéneo, mediante separación de variables, obteniendo

<math>\quad
u(x,t) = \sum_{m=0}^\infty B_{m} cos\left(\frac{m\pi}{b} x\right) e^{-\dfrac{m^2\pi^2}{b^2} t},\quad \text{ donde } \quad
B_m = \frac{\int_{-b}^{b} u(x,0) \cos\left(\frac{m\pi}{b} x\right)dx}{\int_{-b}^{b} \cos^2\left(\frac{m\pi}{b}x\right)dx}.
</math>

Un último detalle antes de proceder a los errores es la obtención simbólica de la solución no acotada.

<math>\quad
\begin{align}
\begin{cases} u_{t} - u_{xx} = 0 \quad &\forall t>0, x\in \mathbb{R},\\
u(x,0) = e^{-x^2} \quad &\forall x\in \mathbb{R}
\end{cases}
\end{align}
\quad
\Longrightarrow
\quad
u(x,t) = \frac{1}{\sqrt{4\pi t}}e^{\frac{-x^2}{4t}} * e^{-x^2} = \frac{e^{\frac{-x^2}{4t+1}}}{\sqrt{4t+1}}.
</math>

Queremos comprobar cuán bien aproximan las soluciones del caso acotado Dirichlet y Neumann a la obtenida por la solución fundamental. Para ello, escojamos un tiempo máximo <math>t_{max} </math> hasta el cuál tomaremos el error en norma <math> L^{\infty} </math>. Queremos estudiar:

* El efecto de aumentar el rango de definición con <math> a,b </math>.
* El efecto de aumentar el número de términos en las series de los sistemas acotados.

=== Análisis de los errores ===
[[Archivo:GrupoADMR Calor Errores1.jpg|550px|thumb|right| Representación gráfica de las soluciones del problema acotado (<math>a = b = 20</math>) y no acotado truncando las series de Fourier en su término <math>100</math>, para un tiempo máximo de <math>1</math>; y de sus errores. Nótese el orden de magnitud de los errores -pequeño-, los picos generados por los errores numéricos dados por la fórmula del trapecio y la suavidad en la frontera dada por las condiciones en esta.]]

Consideramos los parámetros, donde <math> n </math> es el término de la serie de Fourier donde truncamos:

*<math> a = b = 20, </math>
*<math> t_{max} = 1, </math>
*<math> n = 100. </math>

Con <math> t_{max} </math> fijo podemos observar que al aumentar <math>a, b</math> y <math>n</math> el error entre las soluciones acotadas y no acotada disminuye. Nótese que se estanca a partir de cierto <math>n</math>. De hecho, el error llega a aumentar a partir de ciertos <math> a,b </math> si truncamos demasiado pronto o demasiado tarde por la aproximación del trapecio para particiones del espacio demasiado bastas -no observado por pericia en código final-.

{|style="margin: 0 auto;"
| [[Archivo:GrupoADMR Calor Errores MaxD.jpeg|250px|thumb|upright|Con <math>t_{max}=1</math>, vamos obteniendo errores máximos en escala logarítmica con diferentes longitudes y truncamientos entre la solución del problema no acotado y el acotado Dirichlet. Los errores son bastante pequeños si se toman suficientes términos de la serie y una discretización espacial adecuada (en otro caso, aunque no esté representado en esta versión final, los errores divergen)]]
| [[Archivo:GrupoADMR Calor Errores MaxNe.jpeg|250px|thumb|upright|Análogo a la figura adyacente. Nótese que ambas son buenas aproximaciones con suficiente longitud y términos de la serie.]]
|}

Al ampliar el tiempo con la <math>n</math> y fijada <math>a=b=15</math>, al hacer más grande <math>t_{max}</math> los errores máximos no aumentan. En cambio, si tomamos el máximo error para diferentes tiempos vemos con suficiente precisión que no hay diferencia pero que para truncamientos precoces Neumann tiene ligeramente menor error.

{|style="margin: 0 auto;"
| [[Archivo:GrupoADMR Calor Errores DNo.jpeg|250px|thumb|upright|Con dominio fijo, el error en escala logarítmica entre solución no acotada y sistema Dirichlet decae incluso con el paso del tiempo. Esto encaja con que el estado estacionario de Dirichlet será la constante de temperatura nula, al igual que la convolución con la solución fundamental.]]
| [[Archivo:GrupoADMR Calor Errores NoNe.jpeg|250px|thumb|upright|Con condiciones análogas en frontera Neumann, el resultado es igualmente excelente para aproximar el caso no acotado. Aunque en el caso Neumann la solución estacionaria no sea nula, sino un reparto equitativo del calor inicial, cuando el rango es suficientemente amplio, la aproximación será adecuada.]]
| [[Archivo:GrupoADMR Calor Errores DNe.jpeg|250px|thumb|upright|Podemos ver además que ambos sistemas se asemejan cada vez más, por lo que ambas aproximaciones son igualmente útiles para el caso no acotado.]]
|}

=Coeficiente de difusión=

[[Archivo:Fundamental.gif|300px|thumb|right|Solución fundamental de la ecuación del calor unidimensional para distintos valores del coeficiente de difusión. La representación gráfica parte del tiempo inicial <math/> 0.01, <math/> esperando una delta de Dirac en el instante inicial. Nótese que en este tiempo, a mayor coeficiente de difusión más se ha difundido el calor en el sistema, lo que se nota especialmente en el valor máximo observado y su disminución cuanto más crece el coeficiente.]]

Interpretemos el valor de <math/> D <math/> aprovechando el montaje físico anterior. Numéricamente vemos que cuanto mayor es este coeficiente más rápido se difunde el calor, evidenciado en la caída de la solución en tiempos cada vez más cortos.

[[Archivo:Convoluciones_teorica.gif|500px|thumb|right|Soluciones de la ecuación del calor unidimensional para distintos valores del coeficiente de difusión.]]

En el caso unidimensional acotado con solución

<math>
u(x,t) = \sum_{n=1}^{\infty} A_n e^{-D (n\pi)^2 t} \sin(n\pi x).
</math>

Podemos ver que un mayor coeficiente <math/> D <math/> acelera el decaimiento de la exponencial por lo que tendrá un efecto similar justificando su nombre.

= Códigos =

<source lang="matlab">
%%% Código generado junto con Chat GPT para la representación de la solución fundamental de la ecuación del calor unidimensional para distintos valores del coeficiente de difusión

clc
clear all
close all

% Vector de valores de D (coeficiente de difusión)
DD = [1,2,3,4,5,6,7,8,9,10];

% Nombre del archivo GIF de salida
output_gif = 'Fundamental.gif';

% Crear la figura
figura = figure(1);
grid on;
view(3);

% Bucle para iterar sobre cada valor de D
for i = 1:length(DD)
D = DD(i);

% Función fundamental de la ecuación del calor
phi_D = @(x, t) exp(-x.^2./(4*D*t))./sqrt(4*D*pi*t);

% Definición de la malla
xx = -1:0.01:1; % Rango de x
tt = 10^(-2):0.01:1; % Rango de t
[X, T] = meshgrid(xx, tt); % Malla de coordenadas
PHI_D = phi_D(X, T); % Evaluación de la función

% Graficar la superficie
clf;
surf(X, T, PHI_D);
shading flat
shading interp
colormap('jet')
xlabel('x');
ylabel('t');
zlabel('$\Phi_D(x,t)$', 'Interpreter', 'latex');
title(['Solución fundamental de la ecuación del calor, cuando D=', num2str(D)])
axis([-1, 1, 0, 1, 0, 3])
drawnow;

% Captura de la imagen para el GIF
frame = getframe(figura);
img = frame2im(frame);
[imind, cm] = rgb2ind(img, 256);

% Escribir la imagen en el archivo GIF
if i == 1
imwrite(imind, cm, output_gif, 'gif', 'LoopCount', Inf, 'DelayTime', 0.7);
else
imwrite(imind, cm, output_gif, 'gif', 'WriteMode', 'append', 'DelayTime', 0.7);
end

end
</source>

[[Archivo:Convoluciones_teorica.gif|500px|thumb|right|Soluciones de la ecuación del calor unidimensional para distintos valores del coeficiente de difusión.]]

<source lang="matlab">
%%% Código generado junto con Chat GPT para la representación de la soluciones de la ecuación del calor unidimensional para distintos valores del coeficiente de difusión

clc
clear all
close all

% Vector de valores de D (coeficiente de difusión)
DD = [1,10,20,30,40,50,60,70,80,90,100];

% Función gaussiana
gaus = @(x) (exp(-x.^2));

% Nombre del archivo GIF de salida
output_gif = 'Convoluciones_teorica.gif';

% Crear la figura
figura = figure(1);
grid on;
view(3);

% Bucle para iterar sobre cada valor de D
for h = 1:length(DD)
D = DD(h);

% Funciones fundamentales y convolución
phi_D = @(x, t) exp(-x.^2./(4*D*t))./sqrt(4*D*pi*t);
convo = @(x,t) (exp(-x.^2./(4*D.*t+1))/sqrt(4*D.*t+1));

% Definición de la malla
xx = -10:0.1:10; % Rango de x
tt = 10^(-2):0.01:1; % Rango de t

% Inicialización de la matriz de convolución
conv = zeros(length(xx), length(tt));

% Cálculo de la convolución
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
val_conv(i,j) = convo(auxx, auxt);
end
end

% Graficar la superficie
clf;
hold on;
view(3);
surf(xx, tt, val_conv.');
xlabel('x');
ylabel('t');
zlabel('$\Phi_D(x,t)*N(x)$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap('jet');
title(['Soluciones en dominio no acotado, cuando D=', num2str(D)]);
axis([-10, 10, 0, 1, 0, 1]);
hold off;
drawnow;

% Captura de la imagen para el GIF
frame = getframe(figura);
img = frame2im(frame);
[imind, cm] = rgb2ind(img, 256);

% Escribir la imagen en el archivo GIF
if h == 1
imwrite(imind, cm, output_gif, 'gif', 'LoopCount', Inf, 'DelayTime', 0.7);
else
imwrite(imind, cm, output_gif, 'gif', 'WriteMode', 'append', 'DelayTime', 0.7);
end

end
</source>

<source lang="matlab">
%%% Código generado junto con Chat GPT para la representación de los errores entre las soluciones

function [Diff12,Difft1,Difft2] = errores2(n_max, a, b,t_max)

% Función como condición inicial: una función Gaussiana
gaus = @(x) (exp(-x.^2));
u0 = @(x) exp(-x.^2);
D = 1; % Coeficiente de difusión

% Funciones fundamentales y convolución para el modelo de difusión
phi_D = @(x, t) exp(-x.^2./(4*D*t))./sqrt(4*D*pi*t);
convo = @(x,t) (exp(-x.^2./(4*D.*t+1))/sqrt(4*D.*t+1));

% Definición de la malla (espacio y tiempo)
xx = -a:0.01:a;
tt = 10^(-2):0.01:t_max;
[X, T] = meshgrid(xx, tt);

% Inicialización de la matriz de convolución
conv = zeros(length(xx), length(tt));

% Cálculo de la convolución
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
val_conv(j,i) = convo(auxx, auxt);
end
end

% Inicialización de la matriz de resultados para la solución Dirichlet
res_dir = zeros(length(xx), length(tt));

% Cálculo de la solución directa
for n=1:n_max
A(n) = trapz(xx, u0(xx) .* cos((pi/a)*(n - 1/2) * xx)) /...
trapz(xx, cos((pi/a)*(n - 1/2) * xx).^2);
end

% Evaluación de la solución directa en el dominio espacio-temporal
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
u = 0;

% Sumar los términos de la serie de Fourier
for n=1:n_max
u = u + A(n) .* cos((pi/a)*(n - 1/2) * auxx) .*...
exp(-((pi^2 / a^2) * (n - 1/2)^2) .* auxt);
end
res_dir(i,j) = u;
end
end

% Inicialización de la matriz de resultados para la solución en frontera Neumann
res_neu = zeros(length(xx), length(tt));

% Cálculo de la solución en frontera Neumann
for n=1:n_max
B(n) = trapz(xx, u0(xx) .* cos((n*pi/b) * xx)) /...
trapz(xx, cos((n*pi/b) * xx).^2);
end

% Evaluación de la solución en frontera Neumann en el dominio espacio-temporal
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
u = trapz(xx, u0(xx))./(2.*b);

% Sumar los términos de la serie de Fourier
for n=1:n_max
u = u + B(n) .* cos((n*pi/b) * auxx) .* exp(-((n^2 * pi^2 / b^2)) .* auxt);
end
res_neu(i,j) = u;
end
end

% Transponer los resultados para compararlos con la convolución
U1 = res_dir.';
U2 = res_neu.';

% Calcular las diferencias entre las soluciones y la convolución
Diff12 = abs(U1 - U2);
Difft1 = abs(val_conv - U1);
Difft2 = abs(val_conv - U2);

end
</source>

<source lang="matlab">
%%% Código generado junto con Chat GPT para la representación de las soluciones de la ecuación del calor unidimensional junto con la comparación de sus errores

clc
clear all
close all

% Función como condicion inicial
% Se define la función inicial gausiana y la función u0 que se utilizarán más adelante para los cálculos.
gaus = @(x) (exp(-x.^2));
u0 = @(x) exp(-x.^2);
D = 1;

% Funciones fundamentales y convolución
% Se definen las funciones que describen la solución fundamental de la ecuación y la convolución utilizada.
phi_D = @(x, t) exp(-x.^2./(4*D*t))./sqrt(4*D*pi*t);
convo = @(x,t) (exp(-x.^2./(4*D.*t+1))/sqrt(4*D.*t+1));

% Definición de la malla
% Se definen los parámetros de la malla para las variables espaciales (xx) y temporales (tt), y se crea la malla 2D.
a = 100;
b = 100;
tmax = 1;
n_max = 1000;

xx = linspace(-a, a, 1000);
tt = linspace(10^(-2), tmax, 1000);
[X, T] = meshgrid(xx, tt);

% Inicialización de la matriz de convolución
% Se crea una matriz vacía para almacenar los valores de la convolución.
conv = zeros(length(xx), length(tt));

% Cálculo de la convolución
% Se calcula la convolución para cada punto en la malla de la variable temporal tt y espacial xx.
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
val_conv(i,j) = convo(auxx, auxt);
end
end

% Resolución de la ecuación con condiciones de frontera Dirichlet
% Se calculan los coeficientes A para la expansión en series de Fourier y se calcula la solución directa en cada paso temporal.
res_dir = zeros(length(xx), length(tt));

for n=1:n_max
A(n) = trapz(xx, u0(xx) .* cos((pi/a)*(n - 1/2) * xx)) /...
trapz(xx, cos((pi/a)*(n - 1/2) * xx).^2);
end

for j = 1:length(tt)
auxt = tt(j);

for i = 1:length(xx)
auxx = xx(i);
u = 0;

for n=1:n_max
u = u + A(n) .* cos((pi/a)*(n - 1/2) * auxx) .*...
exp(-((pi^2 / a^2) * (n - 1/2)^2) .* auxt);

end
res_dir(j,i) = u;
end
end

% Resolución de la ecuación con condiciones de frontera Neumann
% Se calculan los coeficientes B para la expansión en series de Fourier con las condiciones de frontera Neumann y se calcula la solución en cada paso temporal.
res_neu = zeros(length(xx), length(tt));

for n=1:n_max
B(n) = trapz(xx, u0(xx) .* cos((n*pi/b) * xx)) /...
trapz(xx, cos((n*pi/b) * xx).^2);
end

for j = 1:length(tt)
auxt = tt(j);

for i = 1:length(xx)
auxx = xx(i);
u = trapz(xx, u0(xx))./(2.*b);

for n=1:n_max
u = u + B(n) .* cos((n*pi/b) * auxx) .* exp(-((n^2 * pi^2 / b^2)) .* auxt);

end
res_neu(j,i) = u;
end
end

% Comparación entre las soluciones obtenidas
% Se calcula la diferencia entre las soluciones directa y Neumann, así como entre las soluciones y la convolución.
U1 = res_dir;
U2 = res_neu;

Diff12 = abs(U1 - U2);
Difft1 = abs(val_conv.' - U1);
Difft2 = abs(val_conv.' - U2);

% Visualización de los resultados
% Se muestran las superficies 3D de las soluciones y sus diferencias en gráficos para comparar los resultados obtenidos.
figure(1);

subplot(2,4,2);
view(3);
surf(X, T, U1);
xlabel('x');
ylabel('t');
zlabel('$u_1(x,t)$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca,'jet');
axis([-a,a,0,tmax,0,1])
title('$u_1(x,t)$', 'Interpreter', 'latex');

subplot(2,4,3);
view(3);
surf(X, T, U2);
xlabel('x');
ylabel('t');
zlabel('$u_2(x,t)$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca, 'jet');
axis([-b,b,0,tmax,0,1])
title('$u_2(x,t)$', 'Interpreter', 'latex');

subplot(2,4,4);
view(3);
surf(X, T, Diff12);
xlabel('x');
ylabel('t');
zlabel('$|u_1(x,t)-u_2(x,t)|$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca, 'hot');
title('$|u_1(x,t)-u_2(x,t)|$', 'Interpreter', 'latex');

subplot(2,4,5)
view(3);
surf(X, T, val_conv.');
xlabel('x');
ylabel('t');
zlabel('$\Phi_D(x,t)*N(x)$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca, 'jet');
axis([-a,a,0,tmax,0,1])
title('$\Phi_D(x,t)*N(x)$', 'Interpreter', 'latex');

subplot(2,4,6);
view(3);
surf(X, T, Difft1);
xlabel('x');
ylabel('t');
zlabel('$|\Phi_D(x,t)*N(x)-u_1(x,t)|$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca, 'hot');
title('$|\Phi_D(x,t)*N(x)-u_1(x,t)|$', 'Interpreter', 'latex');

subplot(2,4,7);
view(3);
surf(X, T, Difft2);
xlabel('x');
ylabel('t');
zlabel('$|\Phi_D(x,t)*N(x)-u_2(x,t)|$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca, 'hot');
title('$|\Phi_D(x,t)*N(x)-u_2(x,t)|$', 'Interpreter', 'latex');
</source>

= Referencias =

* Sandro Salsa, Gianmaria Verzini. Partial Differential Equations in Action From Modelling to Theory ([https://link.springer.com/book/10.1007/978-3-031-21853-8 [1<nowiki>]</nowiki>]).

[[Categoría:EDP]]
[[Categoría:EDP24/25]]

Ecuación del calor (ADMR)

2025-03-19T07:37:15Z

Daniel Rodríguez Calderón: /* Coeficiente de difusión */

{{ TrabajoED | Ecuación del calor (Grupo ADMR). | [[:Categoría:EDP|EDP]]|[[:Categoría:EDP24/25|2024-25]] |
*Ángel De Lucas Miranda,
*Daniel Rodríguez Calderón,
*Marcos Cabellos Hernández,
*Rafael Pascual Ortega.}}

=Introducción=

El desarrollo de la ecuación del calor se atribuye principalmente al matemático y físico francés Joseph Fourier, quien en 1822 publicó su obra ''Théorie analytique de la chaleur''. En este trabajo formuló una ecuación diferencial en derivadas parciales que describe la propagación del calor por difusión en un medio continuo. En esta ocasión, profundicemos en dos aspectos de este modelo.

Desde un punto de vista físico, estas ecuaciones permiten la modelización de fenómenos como el calentamiento de materiales o la difusión de solutos en fluidos estacionarios. Sin embargo, la influencia de la física es más profunda: de hecho, el coeficiente de difusión lleva su nombre por la interpretación física de su acción en el sistema. ¿Cuál es este efecto, y por qué se lo denota de esta forma?

Planteemos también el paso entre problemas de dominio acotado y los de no acotado. ¿Podemos aproximar las soluciones del sistema no acotado mediante las del acotado? Y de ser así, ¿qué condiciones de frontera en una dimensión funcionan mejor?

=Solución acotada vs Solución no acotada=
Demos por conocidos el uso de separación de variables para resolver problemas de autofunciones en caso acotado y el manejo de la solución fundamental junto con la convolución para el sistema no acotado. Para comparar todas las soluciones correspondientes, consideremos sistema <math> 1-</math>dimensional, con una gaussiana <math> e^{-x^2} </math> como condición inicial. Justificamos físicamente como el uso de un láser para calentar un punto de una tira fina de metal. Considerando que los fotones puedan desviarse como una distribución normal por el teorema central del límite, es un montaje experimental factible y regular. ¿Y qué condiciones de frontera?

Si consideramos de Dirichlet, a una tira de longitud <math> 2a </math> le colocaremos en los extremos un material que se mantenga siempre a la misma temperatura, pongamos que a una temperatura nula. El sistema resulta

<math>\quad
\begin{align}
\begin{cases} u_{t} - u_{xx} = 0 \quad &\forall t > 0, x\in [-a,a],\\
u(-a,t) = u(a,t) = 0 \quad &\forall t > 0, \\
u(x,0) = e^{-x^2} \quad &\forall x\in [-a,a].
\end{cases}
\end{align}
</math>

La solución al resolver el problema de autofunciones viene dada por

<math>\quad
u(x,t) = \sum_{n=1}^\infty A_{n} \cos\left(\frac{\pi}{a}(n - \frac{1}{2}) x\right) e^{-\dfrac{\pi^2}{a^2}\left(n - \frac{1}{2}\right)^2 t},
\quad \text{ donde } \quad
A_n = \frac{\int_{-a}^{a} u(x,0) \cos\left(\frac{\pi}{a}(n - \frac{1}{2}) x\right)dx}{\int_{-a}^{a} \cos^2\left(\frac{\pi}{a}(n - \frac{1}{2}) x\right)dx}.
</math>

Consideremos ahora condiciones de frontera Neumann, con un flujo nulo que físicamente viene de poner extremos aislantes. Para cierto <math> b > 0,</math> y no necesariamente igual que con <math> a </math>. Consideraremos la ecuación del calor

<math>\quad
\begin{align}
\begin{cases} u_{t} - u_{xx} = 0 \quad &\forall t > 0 , x\in [-b,b],\\
u_x(-b,t) = u_x(b,t) = 0 \quad &\forall t > 0, \\
u(x,0) = e^{-x^2} \quad &\forall x\in[-b,b].
\end{cases}
\end{align}
</math>

De nuevo podemos resolver el sistema, por ser homogéneo, mediante separación de variables, obteniendo

<math>\quad
u(x,t) = \sum_{m=0}^\infty B_{m} cos\left(\frac{m\pi}{b} x\right) e^{-\dfrac{m^2\pi^2}{b^2} t},\quad \text{ donde } \quad
B_m = \frac{\int_{-b}^{b} u(x,0) \cos\left(\frac{m\pi}{b} x\right)dx}{\int_{-b}^{b} \cos^2\left(\frac{m\pi}{b}x\right)dx}.
</math>

Un último detalle antes de proceder a los errores es la obtención simbólica de la solución no acotada.

<math>\quad
\begin{align}
\begin{cases} u_{t} - u_{xx} = 0 \quad &\forall t>0, x\in \mathbb{R},\\
u(x,0) = e^{-x^2} \quad &\forall x\in \mathbb{R}
\end{cases}
\end{align}
\quad
\Longrightarrow
\quad
u(x,t) = \frac{1}{\sqrt{4\pi t}}e^{\frac{-x^2}{4t}} * e^{-x^2} = \frac{e^{\frac{-x^2}{4t+1}}}{\sqrt{4t+1}}.
</math>

Queremos comprobar cuán bien aproximan las soluciones del caso acotado Dirichlet y Neumann a la obtenida por la solución fundamental. Para ello, escojamos un tiempo máximo <math>t_{max} </math> hasta el cuál tomaremos el error en norma <math> L^{\infty} </math>. Queremos estudiar:

* El efecto de aumentar el rango de definición con <math> a,b </math>.
* El efecto de aumentar el número de términos en las series de los sistemas acotados.

=== Análisis de los errores ===
[[Archivo:GrupoADMR Calor Errores1.jpg|550px|thumb|right| Representación gráfica de las soluciones del problema acotado (<math>a = b = 20</math>) y no acotado truncando las series de Fourier en su término <math>100</math>, para un tiempo máximo de <math>1</math>; y de sus errores. Nótese el orden de magnitud de los errores -pequeño-, los picos generados por los errores numéricos dados por la fórmula del trapecio y la suavidad en la frontera dada por las condiciones en esta.]]

Consideramos los parámetros, donde <math> n </math> es el término de la serie de Fourier donde truncamos:

*<math> a = b = 20, </math>
*<math> t_{max} = 1, </math>
*<math> n = 100. </math>

Con <math> t_{max} </math> fijo podemos observar que al aumentar <math>a, b</math> y <math>n</math> el error entre las soluciones acotadas y no acotada disminuye. Nótese que se estanca a partir de cierto <math>n</math>. De hecho, el error llega a aumentar a partir de ciertos <math> a,b </math> si truncamos demasiado pronto o demasiado tarde por la aproximación del trapecio para particiones del espacio demasiado bastas -no observado por pericia en código final-.

{|style="margin: 0 auto;"
| [[Archivo:GrupoADMR Calor Errores MaxD.jpeg|250px|thumb|upright|Con <math>t_{max}=1</math>, vamos obteniendo errores máximos en escala logarítmica con diferentes longitudes y truncamientos entre la solución del problema no acotado y el acotado Dirichlet. Los errores son bastante pequeños si se toman suficientes términos de la serie y una discretización espacial adecuada (en otro caso, aunque no esté representado en esta versión final, los errores divergen)]]
| [[Archivo:GrupoADMR Calor Errores MaxNe.jpeg|250px|thumb|upright|Análogo a la figura adyacente. Nótese que ambas son buenas aproximaciones con suficiente longitud y términos de la serie.]]
|}

Al ampliar el tiempo con la <math>n</math> y fijada <math>a=b=15</math>, al hacer más grande <math>t_{max}</math> los errores máximos no aumentan. En cambio, si tomamos el máximo error para diferentes tiempos vemos con suficiente precisión que no hay diferencia pero que para truncamientos precoces Neumann tiene ligeramente menor error.

{|style="margin: 0 auto;"
| [[Archivo:GrupoADMR Calor Errores DNo.jpeg|250px|thumb|upright|Con dominio fijo, el error en escala logarítmica entre solución no acotada y sistema Dirichlet decae incluso con el paso del tiempo. Esto encaja con que el estado estacionario de Dirichlet será la constante de temperatura nula, al igual que la convolución con la solución fundamental.]]
| [[Archivo:GrupoADMR Calor Errores NoNe.jpeg|250px|thumb|upright|Con condiciones análogas en frontera Neumann, el resultado es igualmente excelente para aproximar el caso no acotado. Aunque en el caso Neumann la solución estacionaria no sea nula, sino un reparto equitativo del calor inicial, cuando el rango es suficientemente amplio, la aproximación será adecuada.]]
| [[Archivo:GrupoADMR Calor Errores DNe.jpeg|250px|thumb|upright|Podemos ver además que ambos sistemas se asemejan cada vez más, por lo que ambas aproximaciones son igualmente útiles para el caso no acotado.]]
|}

=Coeficiente de difusión=

[[Archivo:Fundamental.gif|300px|thumb|right|Solución fundamental de la ecuación del calor unidimensional para distintos valores del coeficiente de difusión. La representación gráfica parte del tiempo inicial <math/> 0.01, <math/> esperando una delta de Dirac en el instante inicial. Nótese que en este tiempo, a mayor coeficiente de difusión más se ha difundido el calor en el sistema, lo que se nota especialmente en el valor máximo observado y su disminución cuanto más crece el coeficiente.]]

Interpretemos el valor de <math/> D <math/> aprovechando el montaje físico anterior. Numéricamente vemos que cuanto mayor es este coeficiente más rápido se difunde el calor, evidenciado en la caída de la solución en tiempos cada vez más cortos.

[[Archivo:Convoluciones_teorica.gif|500px|thumb|right|Soluciones de la ecuación del calor unidimensional para distintos valores del coeficiente de difusión.]]

En el caso unidimensional acotado con solución

<math>
u(x,t) = \sum_{n=1}^{\infty} A_n e^{-D (n\pi)^2 t} \sin(n\pi x),
</math>

Podemos ver que un mayor coeficiente <math/> D <math/> acelera el decaimiento de la exponencial por lo que tendrá un efecto similar justificando su nombre.

= Códigos =

<source lang="matlab">
%%% Código generado junto con Chat GPT para la representación de la solución fundamental de la ecuación del calor unidimensional para distintos valores del coeficiente de difusión

clc
clear all
close all

% Vector de valores de D (coeficiente de difusión)
DD = [1,2,3,4,5,6,7,8,9,10];

% Nombre del archivo GIF de salida
output_gif = 'Fundamental.gif';

% Crear la figura
figura = figure(1);
grid on;
view(3);

% Bucle para iterar sobre cada valor de D
for i = 1:length(DD)
D = DD(i);

% Función fundamental de la ecuación del calor
phi_D = @(x, t) exp(-x.^2./(4*D*t))./sqrt(4*D*pi*t);

% Definición de la malla
xx = -1:0.01:1; % Rango de x
tt = 10^(-2):0.01:1; % Rango de t
[X, T] = meshgrid(xx, tt); % Malla de coordenadas
PHI_D = phi_D(X, T); % Evaluación de la función

% Graficar la superficie
clf;
surf(X, T, PHI_D);
shading flat
shading interp
colormap('jet')
xlabel('x');
ylabel('t');
zlabel('$\Phi_D(x,t)$', 'Interpreter', 'latex');
title(['Solución fundamental de la ecuación del calor, cuando D=', num2str(D)])
axis([-1, 1, 0, 1, 0, 3])
drawnow;

% Captura de la imagen para el GIF
frame = getframe(figura);
img = frame2im(frame);
[imind, cm] = rgb2ind(img, 256);

% Escribir la imagen en el archivo GIF
if i == 1
imwrite(imind, cm, output_gif, 'gif', 'LoopCount', Inf, 'DelayTime', 0.7);
else
imwrite(imind, cm, output_gif, 'gif', 'WriteMode', 'append', 'DelayTime', 0.7);
end

end
</source>

[[Archivo:Convoluciones_teorica.gif|500px|thumb|right|Soluciones de la ecuación del calor unidimensional para distintos valores del coeficiente de difusión.]]

<source lang="matlab">
%%% Código generado junto con Chat GPT para la representación de la soluciones de la ecuación del calor unidimensional para distintos valores del coeficiente de difusión

clc
clear all
close all

% Vector de valores de D (coeficiente de difusión)
DD = [1,10,20,30,40,50,60,70,80,90,100];

% Función gaussiana
gaus = @(x) (exp(-x.^2));

% Nombre del archivo GIF de salida
output_gif = 'Convoluciones_teorica.gif';

% Crear la figura
figura = figure(1);
grid on;
view(3);

% Bucle para iterar sobre cada valor de D
for h = 1:length(DD)
D = DD(h);

% Funciones fundamentales y convolución
phi_D = @(x, t) exp(-x.^2./(4*D*t))./sqrt(4*D*pi*t);
convo = @(x,t) (exp(-x.^2./(4*D.*t+1))/sqrt(4*D.*t+1));

% Definición de la malla
xx = -10:0.1:10; % Rango de x
tt = 10^(-2):0.01:1; % Rango de t

% Inicialización de la matriz de convolución
conv = zeros(length(xx), length(tt));

% Cálculo de la convolución
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
val_conv(i,j) = convo(auxx, auxt);
end
end

% Graficar la superficie
clf;
hold on;
view(3);
surf(xx, tt, val_conv.');
xlabel('x');
ylabel('t');
zlabel('$\Phi_D(x,t)*N(x)$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap('jet');
title(['Soluciones en dominio no acotado, cuando D=', num2str(D)]);
axis([-10, 10, 0, 1, 0, 1]);
hold off;
drawnow;

% Captura de la imagen para el GIF
frame = getframe(figura);
img = frame2im(frame);
[imind, cm] = rgb2ind(img, 256);

% Escribir la imagen en el archivo GIF
if h == 1
imwrite(imind, cm, output_gif, 'gif', 'LoopCount', Inf, 'DelayTime', 0.7);
else
imwrite(imind, cm, output_gif, 'gif', 'WriteMode', 'append', 'DelayTime', 0.7);
end

end
</source>

<source lang="matlab">
%%% Código generado junto con Chat GPT para la representación de los errores entre las soluciones

function [Diff12,Difft1,Difft2] = errores2(n_max, a, b,t_max)

% Función como condición inicial: una función Gaussiana
gaus = @(x) (exp(-x.^2));
u0 = @(x) exp(-x.^2);
D = 1; % Coeficiente de difusión

% Funciones fundamentales y convolución para el modelo de difusión
phi_D = @(x, t) exp(-x.^2./(4*D*t))./sqrt(4*D*pi*t);
convo = @(x,t) (exp(-x.^2./(4*D.*t+1))/sqrt(4*D.*t+1));

% Definición de la malla (espacio y tiempo)
xx = -a:0.01:a;
tt = 10^(-2):0.01:t_max;
[X, T] = meshgrid(xx, tt);

% Inicialización de la matriz de convolución
conv = zeros(length(xx), length(tt));

% Cálculo de la convolución
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
val_conv(j,i) = convo(auxx, auxt);
end
end

% Inicialización de la matriz de resultados para la solución Dirichlet
res_dir = zeros(length(xx), length(tt));

% Cálculo de la solución directa
for n=1:n_max
A(n) = trapz(xx, u0(xx) .* cos((pi/a)*(n - 1/2) * xx)) /...
trapz(xx, cos((pi/a)*(n - 1/2) * xx).^2);
end

% Evaluación de la solución directa en el dominio espacio-temporal
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
u = 0;

% Sumar los términos de la serie de Fourier
for n=1:n_max
u = u + A(n) .* cos((pi/a)*(n - 1/2) * auxx) .*...
exp(-((pi^2 / a^2) * (n - 1/2)^2) .* auxt);
end
res_dir(i,j) = u;
end
end

% Inicialización de la matriz de resultados para la solución en frontera Neumann
res_neu = zeros(length(xx), length(tt));

% Cálculo de la solución en frontera Neumann
for n=1:n_max
B(n) = trapz(xx, u0(xx) .* cos((n*pi/b) * xx)) /...
trapz(xx, cos((n*pi/b) * xx).^2);
end

% Evaluación de la solución en frontera Neumann en el dominio espacio-temporal
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
u = trapz(xx, u0(xx))./(2.*b);

% Sumar los términos de la serie de Fourier
for n=1:n_max
u = u + B(n) .* cos((n*pi/b) * auxx) .* exp(-((n^2 * pi^2 / b^2)) .* auxt);
end
res_neu(i,j) = u;
end
end

% Transponer los resultados para compararlos con la convolución
U1 = res_dir.';
U2 = res_neu.';

% Calcular las diferencias entre las soluciones y la convolución
Diff12 = abs(U1 - U2);
Difft1 = abs(val_conv - U1);
Difft2 = abs(val_conv - U2);

end
</source>

<source lang="matlab">
%%% Código generado junto con Chat GPT para la representación de las soluciones de la ecuación del calor unidimensional junto con la comparación de sus errores

clc
clear all
close all

% Función como condicion inicial
% Se define la función inicial gausiana y la función u0 que se utilizarán más adelante para los cálculos.
gaus = @(x) (exp(-x.^2));
u0 = @(x) exp(-x.^2);
D = 1;

% Funciones fundamentales y convolución
% Se definen las funciones que describen la solución fundamental de la ecuación y la convolución utilizada.
phi_D = @(x, t) exp(-x.^2./(4*D*t))./sqrt(4*D*pi*t);
convo = @(x,t) (exp(-x.^2./(4*D.*t+1))/sqrt(4*D.*t+1));

% Definición de la malla
% Se definen los parámetros de la malla para las variables espaciales (xx) y temporales (tt), y se crea la malla 2D.
a = 100;
b = 100;
tmax = 1;
n_max = 1000;

xx = linspace(-a, a, 1000);
tt = linspace(10^(-2), tmax, 1000);
[X, T] = meshgrid(xx, tt);

% Inicialización de la matriz de convolución
% Se crea una matriz vacía para almacenar los valores de la convolución.
conv = zeros(length(xx), length(tt));

% Cálculo de la convolución
% Se calcula la convolución para cada punto en la malla de la variable temporal tt y espacial xx.
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
val_conv(i,j) = convo(auxx, auxt);
end
end

% Resolución de la ecuación con condiciones de frontera Dirichlet
% Se calculan los coeficientes A para la expansión en series de Fourier y se calcula la solución directa en cada paso temporal.
res_dir = zeros(length(xx), length(tt));

for n=1:n_max
A(n) = trapz(xx, u0(xx) .* cos((pi/a)*(n - 1/2) * xx)) /...
trapz(xx, cos((pi/a)*(n - 1/2) * xx).^2);
end

for j = 1:length(tt)
auxt = tt(j);

for i = 1:length(xx)
auxx = xx(i);
u = 0;

for n=1:n_max
u = u + A(n) .* cos((pi/a)*(n - 1/2) * auxx) .*...
exp(-((pi^2 / a^2) * (n - 1/2)^2) .* auxt);

end
res_dir(j,i) = u;
end
end

% Resolución de la ecuación con condiciones de frontera Neumann
% Se calculan los coeficientes B para la expansión en series de Fourier con las condiciones de frontera Neumann y se calcula la solución en cada paso temporal.
res_neu = zeros(length(xx), length(tt));

for n=1:n_max
B(n) = trapz(xx, u0(xx) .* cos((n*pi/b) * xx)) /...
trapz(xx, cos((n*pi/b) * xx).^2);
end

for j = 1:length(tt)
auxt = tt(j);

for i = 1:length(xx)
auxx = xx(i);
u = trapz(xx, u0(xx))./(2.*b);

for n=1:n_max
u = u + B(n) .* cos((n*pi/b) * auxx) .* exp(-((n^2 * pi^2 / b^2)) .* auxt);

end
res_neu(j,i) = u;
end
end

% Comparación entre las soluciones obtenidas
% Se calcula la diferencia entre las soluciones directa y Neumann, así como entre las soluciones y la convolución.
U1 = res_dir;
U2 = res_neu;

Diff12 = abs(U1 - U2);
Difft1 = abs(val_conv.' - U1);
Difft2 = abs(val_conv.' - U2);

% Visualización de los resultados
% Se muestran las superficies 3D de las soluciones y sus diferencias en gráficos para comparar los resultados obtenidos.
figure(1);

subplot(2,4,2);
view(3);
surf(X, T, U1);
xlabel('x');
ylabel('t');
zlabel('$u_1(x,t)$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca,'jet');
axis([-a,a,0,tmax,0,1])
title('$u_1(x,t)$', 'Interpreter', 'latex');

subplot(2,4,3);
view(3);
surf(X, T, U2);
xlabel('x');
ylabel('t');
zlabel('$u_2(x,t)$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca, 'jet');
axis([-b,b,0,tmax,0,1])
title('$u_2(x,t)$', 'Interpreter', 'latex');

subplot(2,4,4);
view(3);
surf(X, T, Diff12);
xlabel('x');
ylabel('t');
zlabel('$|u_1(x,t)-u_2(x,t)|$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca, 'hot');
title('$|u_1(x,t)-u_2(x,t)|$', 'Interpreter', 'latex');

subplot(2,4,5)
view(3);
surf(X, T, val_conv.');
xlabel('x');
ylabel('t');
zlabel('$\Phi_D(x,t)*N(x)$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca, 'jet');
axis([-a,a,0,tmax,0,1])
title('$\Phi_D(x,t)*N(x)$', 'Interpreter', 'latex');

subplot(2,4,6);
view(3);
surf(X, T, Difft1);
xlabel('x');
ylabel('t');
zlabel('$|\Phi_D(x,t)*N(x)-u_1(x,t)|$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca, 'hot');
title('$|\Phi_D(x,t)*N(x)-u_1(x,t)|$', 'Interpreter', 'latex');

subplot(2,4,7);
view(3);
surf(X, T, Difft2);
xlabel('x');
ylabel('t');
zlabel('$|\Phi_D(x,t)*N(x)-u_2(x,t)|$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca, 'hot');
title('$|\Phi_D(x,t)*N(x)-u_2(x,t)|$', 'Interpreter', 'latex');
</source>

= Referencias =

* Sandro Salsa, Gianmaria Verzini. Partial Differential Equations in Action From Modelling to Theory ([https://link.springer.com/book/10.1007/978-3-031-21853-8 [1<nowiki>]</nowiki>]).

[[Categoría:EDP]]
[[Categoría:EDP24/25]]

Ecuación del calor (ADMR)

2025-03-18T23:04:49Z

Daniel Rodríguez Calderón: /* Códigos */

{{ TrabajoED | Ecuación del calor (Grupo ADMR). | [[:Categoría:EDP|EDP]]|[[:Categoría:EDP24/25|2024-25]] |
*Ángel De Lucas Miranda,
*Daniel Rodríguez Calderón,
*Marcos Cabellos Hernández,
*Rafael Pascual Ortega.}}

=Introducción=

El desarrollo de la ecuación del calor se atribuye principalmente al matemático y físico francés Joseph Fourier, quien en 1822 publicó su obra ''Théorie analytique de la chaleur''. En este trabajo formuló una ecuación diferencial en derivadas parciales que describe la propagación del calor por difusión en un medio continuo. En esta ocasión, profundicemos en dos aspectos de este modelo.

Desde un punto de vista físico, estas ecuaciones permiten la modelización de fenómenos como el calentamiento de materiales o la difusión de solutos en fluidos estacionarios. Sin embargo, la influencia de la física es más profunda: de hecho, el coeficiente de difusión lleva su nombre por la interpretación física de su acción en el sistema. ¿Cuál es este efecto, y por qué se lo denota de esta forma?

Planteemos también el paso entre problemas de dominio acotado y los de no acotado. ¿Podemos aproximar las soluciones del sistema no acotado mediante las del acotado? Y de ser así, ¿qué condiciones de frontera en una dimensión funcionan mejor?

=Solución acotada vs Solución no acotada=
Demos por conocidos el uso de separación de variables para resolver problemas de autofunciones en caso acotado y el manejo de la solución fundamental junto con la convolución para el sistema no acotado. Para comparar todas las soluciones correspondientes, consideremos sistema <math> 1-</math>dimensional, con una gaussiana <math> e^{-x^2} </math> como condición inicial. Justificamos físicamente como el uso de un láser para calentar un punto de una tira fina de metal. Considerando que los fotones puedan desviarse como una distribución normal por el teorema central del límite, es un montaje experimental factible y regular. ¿Y qué condiciones de frontera?

Si consideramos de Dirichlet, a una tira de longitud <math> 2a </math> le colocaremos en los extremos un material que se mantenga siempre a la misma temperatura, pongamos que a una temperatura nula. El sistema resulta

<math>\quad
\begin{align}
\begin{cases} u_{t} - u_{xx} = 0 \quad &\forall t > 0, x\in [-a,a],\\
u(-a,t) = u(a,t) = 0 \quad &\forall t > 0, \\
u(x,0) = e^{-x^2} \quad &\forall x\in [-a,a].
\end{cases}
\end{align}
</math>

La solución al resolver el problema de autofunciones viene dada por

<math>\quad
u(x,t) = \sum_{n=1}^\infty A_{n} \cos\left(\frac{\pi}{a}(n - \frac{1}{2}) x\right) e^{-\dfrac{\pi^2}{a^2}\left(n - \frac{1}{2}\right)^2 t},
\quad \text{ donde } \quad
A_n = \frac{\int_{-a}^{a} u(x,0) \cos\left(\frac{\pi}{a}(n - \frac{1}{2}) x\right)dx}{\int_{-a}^{a} \cos^2\left(\frac{\pi}{a}(n - \frac{1}{2}) x\right)dx}.
</math>

Consideremos ahora condiciones de frontera Neumann, con un flujo nulo que físicamente viene de poner extremos aislantes. Para cierto <math> b > 0,</math> y no necesariamente igual que con <math> a </math>. Consideraremos la ecuación del calor

<math>\quad
\begin{align}
\begin{cases} u_{t} - u_{xx} = 0 \quad &\forall t > 0 , x\in [-b,b],\\
u_x(-b,t) = u_x(b,t) = 0 \quad &\forall t > 0, \\
u(x,0) = e^{-x^2} \quad &\forall x\in[-b,b].
\end{cases}
\end{align}
</math>

De nuevo podemos resolver el sistema, por ser homogéneo, mediante separación de variables, obteniendo

<math>\quad
u(x,t) = \sum_{m=0}^\infty B_{m} cos\left(\frac{m\pi}{b} x\right) e^{-\dfrac{m^2\pi^2}{b^2} t},\quad \text{ donde } \quad
B_m = \frac{\int_{-b}^{b} u(x,0) \cos\left(\frac{m\pi}{b} x\right)dx}{\int_{-b}^{b} \cos^2\left(\frac{m\pi}{b}x\right)dx}.
</math>

Un último detalle antes de proceder a los errores es la obtención simbólica de la solución no acotada.

<math>\quad
\begin{align}
\begin{cases} u_{t} - u_{xx} = 0 \quad &\forall t>0, x\in \mathbb{R},\\
u(x,0) = e^{-x^2} \quad &\forall x\in \mathbb{R}
\end{cases}
\end{align}
\quad
\Longrightarrow
\quad
u(x,t) = \frac{1}{\sqrt{4\pi t}}e^{\frac{-x^2}{4t}} * e^{-x^2} = \frac{e^{\frac{-x^2}{4t+1}}}{\sqrt{4t+1}}.
</math>

Queremos comprobar cuán bien aproximan las soluciones del caso acotado Dirichlet y Neumann a la obtenida por la solución fundamental. Para ello, escojamos un tiempo máximo <math>t_{max} </math> hasta el cuál tomaremos el error en norma <math> L^{\infty} </math>. Queremos estudiar:

* El efecto de aumentar el rango de definición con <math> a,b </math>.
* El efecto de aumentar el número de términos en las series de los sistemas acotados.

=== Análisis de los errores ===
[[Archivo:GrupoADMR Calor Errores1.jpg|550px|thumb|right| Representación gráfica de las soluciones del problema acotado (<math>a = b = 20</math>) y no acotado truncando las series de Fourier en su término <math>100</math>, para un tiempo máximo de <math>1</math>; y de sus errores. Nótese el orden de magnitud de los errores -pequeño-, los picos generados por los errores numéricos dados por la fórmula del trapecio y la suavidad en la frontera dada por las condiciones en esta.]]

Consideramos los parámetros, donde <math> n </math> es el término de la serie de Fourier donde truncamos:

*<math> a = b = 20, </math>
*<math> t_{max} = 1, </math>
*<math> n = 100. </math>

Con <math> t_{max} </math> fijo podemos observar que al aumentar <math>a, b</math> y <math>n</math> el error entre las soluciones acotadas y no acotada disminuye. Nótese que se estanca a partir de cierto <math>n</math>. De hecho, el error llega a aumentar a partir de ciertos <math> a,b </math> si truncamos demasiado pronto o demasiado tarde por la aproximación del trapecio para particiones del espacio demasiado bastas -no observado por pericia en código final-.

{|style="margin: 0 auto;"
| [[Archivo:GrupoADMR Calor Errores MaxD.jpeg|250px|thumb|upright|]]
| [[Archivo:GrupoADMR Calor Errores MaxNe.jpeg|250px|thumb|upright|]]
|}

Al ampliar el tiempo con la <math>n</math> y fijada la longitud, al hacer más grande <math>t_{max}</math> los errores máximos no aumentan. En cambio, si tomamos el máximo error para diferentes tiempos vemos con suficiente precisión que no hay diferencia pero que para truncamientos precoces Neumann tiene ligeramente menor error.

{|style="margin: 0 auto;"
| [[Archivo:GrupoADMR Calor Errores DNo.jpeg|250px|thumb|upright|]]
| [[Archivo:GrupoADMR Calor Errores NoNe.jpeg|250px|thumb|upright|]]
| [[Archivo:GrupoADMR Calor Errores DNe.jpeg|250px|thumb|upright|]]
|}

=Coeficiente de difusión=

[[Archivo:Fundamental.gif|300px|thumb|right|Solución fundamental de la ecuación del calor unidimensional para distintos valores del coeficiente de difusión. La representación gráfica parte del tiempo inicial <math/> 0.01, <math/> esperando una delta de Dirac en el instante inicial. Nótese que en este tiempo, a mayor coeficiente de difusión más se ha difundido el calor en el sistema, lo que se nota especialmente en el valor máximo observado y su disminución cuanto más crece el coeficiente.]]

Vamos a interpretar como influye exactamente el coeficiente de difusión en la solución de la ecuación del calor. Imaginemos por un segundo que tenemos varias barras (cada una de un material diferente) y queremos estudiar como se difunde el calor a través de cada una de ellas. Usando Matlab podemos calcular y representar fácilmente la solución del problema para diferentes valores del coeficiente de difusión, que denotaremos por <math/> D <math/>.

## Interpretemos el valor de <math/> D <math/> aprovechando el montaje físico anterior. Numéricamente vemos que cuanto mayor es este coeficiente más rápido se difunde el calor, evidenciado en la caída de la solución en tiempos cada vez más cortos.

[[Archivo:Fundamental.gif|400px|thumb|right|Solución fundamental de la ecuación del calor unidimensional para distintos valores del coeficiente de difusión. La representación gráfica parte del tiempo inicial <math/> 0.01, <math/> esperando una delta de Dirac en el instante inicial. Nótese que en este tiempo, a mayor coeficiente de difusión más se ha difundido el calor en el sistema, lo que se nota especialmente en el valor máximo observado y su disminución cuanto más crece el coeficiente.]]

Como se puede observar en la animación, fijado un tiempo, a medida que crece el valor del coeficiente de difusión el pico de la solución fundamental se va haciendo más cada vez más pequeño, es decir, a mayor coeficiente de difusión, menor es el tiempo que tarda en expandirse el calor. De hecho, podemos observar este hecho directamente de la expresión de la solución fundamental del calor. La solución de la ecuación del calor unidimensional en una barra de longitud $ L = 1 $, con coeficiente de difusión $ D $ y condición inicial dada por una función gaussiana $ f(x) = e^{-x^2} $, se expresa como una serie de Fourier:

## En el caso unidimensional acotado con solución

<math>
u(x,t) = \sum_{n=1}^{\infty} A_n e^{-D (n\pi)^2 t} \sin(n\pi x),
</math>

[Solución que depende de D]- [Viendo la solución del apartado 3 me estoy dando cuenta que también depende de la longitud de la barra, si estuviésemos trabajando con el problema acotado podríamos ver como afectan estos dos parámetros al tiempo de difusión]

Tomando el límite de D tendiendo a <math/> \infty <math/> vemos que u(x,t) tiende a 0, fijado un t determinado. Por tanto, podemos concluir que la solución del calor alcanzará el estado estacionario antes si su coeficiente del calor es mayor.

## Podemos ver que un mayor coeficiente <math/> D <math/> acelera el decaimiento de la exponencial por lo que tendrá un efecto similar justificando su nombre.

= Códigos =

<source lang="matlab">
%%% Código generado junto con Chat GPT para la representación de la solución fundamental de la ecuación del calor unidimensional para distintos valores del coeficiente de difusión

clc
clear all
close all

% Vector de valores de D (coeficiente de difusión)
DD = [1,2,3,4,5,6,7,8,9,10];

% Nombre del archivo GIF de salida
output_gif = 'Fundamental.gif';

% Crear la figura
figura = figure(1);
grid on;
view(3);

% Bucle para iterar sobre cada valor de D
for i = 1:length(DD)
D = DD(i);

% Función fundamental de la ecuación del calor
phi_D = @(x, t) exp(-x.^2./(4*D*t))./sqrt(4*D*pi*t);

% Definición de la malla
xx = -1:0.01:1; % Rango de x
tt = 10^(-2):0.01:1; % Rango de t
[X, T] = meshgrid(xx, tt); % Malla de coordenadas
PHI_D = phi_D(X, T); % Evaluación de la función

% Graficar la superficie
clf;
surf(X, T, PHI_D);
shading flat
shading interp
colormap('jet')
xlabel('x');
ylabel('t');
zlabel('$\Phi_D(x,t)$', 'Interpreter', 'latex');
title(['Solución fundamental de la ecuación del calor, cuando D=', num2str(D)])
axis([-1, 1, 0, 1, 0, 3])
drawnow;

% Captura de la imagen para el GIF
frame = getframe(figura);
img = frame2im(frame);
[imind, cm] = rgb2ind(img, 256);

% Escribir la imagen en el archivo GIF
if i == 1
imwrite(imind, cm, output_gif, 'gif', 'LoopCount', Inf, 'DelayTime', 0.7);
else
imwrite(imind, cm, output_gif, 'gif', 'WriteMode', 'append', 'DelayTime', 0.7);
end

end
</source>

[[Archivo:Convoluciones_teorica.gif|500px|thumb|right|Soluciones de la ecuación del calor unidimensional para distintos valores del coeficiente de difusión.]]

<source lang="matlab">
%%% Código generado junto con Chat GPT para la representación de la soluciones de la ecuación del calor unidimensional para distintos valores del coeficiente de difusión

clc
clear all
close all

% Vector de valores de D (coeficiente de difusión)
DD = [1,10,20,30,40,50,60,70,80,90,100];

% Función gaussiana
gaus = @(x) (exp(-x.^2));

% Nombre del archivo GIF de salida
output_gif = 'Convoluciones_teorica.gif';

% Crear la figura
figura = figure(1);
grid on;
view(3);

% Bucle para iterar sobre cada valor de D
for h = 1:length(DD)
D = DD(h);

% Funciones fundamentales y convolución
phi_D = @(x, t) exp(-x.^2./(4*D*t))./sqrt(4*D*pi*t);
convo = @(x,t) (exp(-x.^2./(4*D.*t+1))/sqrt(4*D.*t+1));

% Definición de la malla
xx = -10:0.1:10; % Rango de x
tt = 10^(-2):0.01:1; % Rango de t

% Inicialización de la matriz de convolución
conv = zeros(length(xx), length(tt));

% Cálculo de la convolución
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
val_conv(i,j) = convo(auxx, auxt);
end
end

% Graficar la superficie
clf;
hold on;
view(3);
surf(xx, tt, val_conv.');
xlabel('x');
ylabel('t');
zlabel('$\Phi_D(x,t)*N(x)$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap('jet');
title(['Soluciones en dominio no acotado, cuando D=', num2str(D)]);
axis([-10, 10, 0, 1, 0, 1]);
hold off;
drawnow;

% Captura de la imagen para el GIF
frame = getframe(figura);
img = frame2im(frame);
[imind, cm] = rgb2ind(img, 256);

% Escribir la imagen en el archivo GIF
if h == 1
imwrite(imind, cm, output_gif, 'gif', 'LoopCount', Inf, 'DelayTime', 0.7);
else
imwrite(imind, cm, output_gif, 'gif', 'WriteMode', 'append', 'DelayTime', 0.7);
end

end
</source>

<source lang="matlab">
%%% Código generado junto con Chat GPT para la representación de los errores entre las soluciones

function [Diff12,Difft1,Difft2] = errores2(n_max, a, b,t_max)

% Función como condición inicial: una función Gaussiana
gaus = @(x) (exp(-x.^2));
u0 = @(x) exp(-x.^2);
D = 1; % Coeficiente de difusión

% Funciones fundamentales y convolución para el modelo de difusión
phi_D = @(x, t) exp(-x.^2./(4*D*t))./sqrt(4*D*pi*t);
convo = @(x,t) (exp(-x.^2./(4*D.*t+1))/sqrt(4*D.*t+1));

% Definición de la malla (espacio y tiempo)
xx = -a:0.01:a;
tt = 10^(-2):0.01:t_max;
[X, T] = meshgrid(xx, tt);

% Inicialización de la matriz de convolución
conv = zeros(length(xx), length(tt));

% Cálculo de la convolución
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
val_conv(j,i) = convo(auxx, auxt);
end
end

% Inicialización de la matriz de resultados para la solución Dirichlet
res_dir = zeros(length(xx), length(tt));

% Cálculo de la solución directa
for n=1:n_max
A(n) = trapz(xx, u0(xx) .* cos((pi/a)*(n - 1/2) * xx)) /...
trapz(xx, cos((pi/a)*(n - 1/2) * xx).^2);
end

% Evaluación de la solución directa en el dominio espacio-temporal
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
u = 0;

% Sumar los términos de la serie de Fourier
for n=1:n_max
u = u + A(n) .* cos((pi/a)*(n - 1/2) * auxx) .*...
exp(-((pi^2 / a^2) * (n - 1/2)^2) .* auxt);
end
res_dir(i,j) = u;
end
end

% Inicialización de la matriz de resultados para la solución en frontera Neumann
res_neu = zeros(length(xx), length(tt));

% Cálculo de la solución en frontera Neumann
for n=1:n_max
B(n) = trapz(xx, u0(xx) .* cos((n*pi/b) * xx)) /...
trapz(xx, cos((n*pi/b) * xx).^2);
end

% Evaluación de la solución en frontera Neumann en el dominio espacio-temporal
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
u = trapz(xx, u0(xx))./(2.*b);

% Sumar los términos de la serie de Fourier
for n=1:n_max
u = u + B(n) .* cos((n*pi/b) * auxx) .* exp(-((n^2 * pi^2 / b^2)) .* auxt);
end
res_neu(i,j) = u;
end
end

% Transponer los resultados para compararlos con la convolución
U1 = res_dir.';
U2 = res_neu.';

% Calcular las diferencias entre las soluciones y la convolución
Diff12 = abs(U1 - U2);
Difft1 = abs(val_conv - U1);
Difft2 = abs(val_conv - U2);

end
</source>

<source lang="matlab">
%%% Código generado junto con Chat GPT para la representación de las soluciones de la ecuación del calor unidimensional junto con la comparación de sus errores

clc
clear all
close all

% Función como condicion inicial
% Se define la función inicial gausiana y la función u0 que se utilizarán más adelante para los cálculos.
gaus = @(x) (exp(-x.^2));
u0 = @(x) exp(-x.^2);
D = 1;

% Funciones fundamentales y convolución
% Se definen las funciones que describen la solución fundamental de la ecuación y la convolución utilizada.
phi_D = @(x, t) exp(-x.^2./(4*D*t))./sqrt(4*D*pi*t);
convo = @(x,t) (exp(-x.^2./(4*D.*t+1))/sqrt(4*D.*t+1));

% Definición de la malla
% Se definen los parámetros de la malla para las variables espaciales (xx) y temporales (tt), y se crea la malla 2D.
a = 100;
b = 100;
tmax = 1;
n_max = 1000;

xx = linspace(-a, a, 1000);
tt = linspace(10^(-2), tmax, 1000);
[X, T] = meshgrid(xx, tt);

% Inicialización de la matriz de convolución
% Se crea una matriz vacía para almacenar los valores de la convolución.
conv = zeros(length(xx), length(tt));

% Cálculo de la convolución
% Se calcula la convolución para cada punto en la malla de la variable temporal tt y espacial xx.
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
val_conv(i,j) = convo(auxx, auxt);
end
end

% Resolución de la ecuación con condiciones de frontera Dirichlet
% Se calculan los coeficientes A para la expansión en series de Fourier y se calcula la solución directa en cada paso temporal.
res_dir = zeros(length(xx), length(tt));

for n=1:n_max
A(n) = trapz(xx, u0(xx) .* cos((pi/a)*(n - 1/2) * xx)) /...
trapz(xx, cos((pi/a)*(n - 1/2) * xx).^2);
end

for j = 1:length(tt)
auxt = tt(j);

for i = 1:length(xx)
auxx = xx(i);
u = 0;

for n=1:n_max
u = u + A(n) .* cos((pi/a)*(n - 1/2) * auxx) .*...
exp(-((pi^2 / a^2) * (n - 1/2)^2) .* auxt);

end
res_dir(j,i) = u;
end
end

% Resolución de la ecuación con condiciones de frontera Neumann
% Se calculan los coeficientes B para la expansión en series de Fourier con las condiciones de frontera Neumann y se calcula la solución en cada paso temporal.
res_neu = zeros(length(xx), length(tt));

for n=1:n_max
B(n) = trapz(xx, u0(xx) .* cos((n*pi/b) * xx)) /...
trapz(xx, cos((n*pi/b) * xx).^2);
end

for j = 1:length(tt)
auxt = tt(j);

for i = 1:length(xx)
auxx = xx(i);
u = trapz(xx, u0(xx))./(2.*b);

for n=1:n_max
u = u + B(n) .* cos((n*pi/b) * auxx) .* exp(-((n^2 * pi^2 / b^2)) .* auxt);

end
res_neu(j,i) = u;
end
end

% Comparación entre las soluciones obtenidas
% Se calcula la diferencia entre las soluciones directa y Neumann, así como entre las soluciones y la convolución.
U1 = res_dir;
U2 = res_neu;

Diff12 = abs(U1 - U2);
Difft1 = abs(val_conv.' - U1);
Difft2 = abs(val_conv.' - U2);

% Visualización de los resultados
% Se muestran las superficies 3D de las soluciones y sus diferencias en gráficos para comparar los resultados obtenidos.
figure(1);

subplot(2,4,2);
view(3);
surf(X, T, U1);
xlabel('x');
ylabel('t');
zlabel('$u_1(x,t)$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca,'jet');
axis([-a,a,0,tmax,0,1])
title('$u_1(x,t)$', 'Interpreter', 'latex');

subplot(2,4,3);
view(3);
surf(X, T, U2);
xlabel('x');
ylabel('t');
zlabel('$u_2(x,t)$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca, 'jet');
axis([-b,b,0,tmax,0,1])
title('$u_2(x,t)$', 'Interpreter', 'latex');

subplot(2,4,4);
view(3);
surf(X, T, Diff12);
xlabel('x');
ylabel('t');
zlabel('$|u_1(x,t)-u_2(x,t)|$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca, 'hot');
title('$|u_1(x,t)-u_2(x,t)|$', 'Interpreter', 'latex');

subplot(2,4,5)
view(3);
surf(X, T, val_conv.');
xlabel('x');
ylabel('t');
zlabel('$\Phi_D(x,t)*N(x)$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca, 'jet');
axis([-a,a,0,tmax,0,1])
title('$\Phi_D(x,t)*N(x)$', 'Interpreter', 'latex');

subplot(2,4,6);
view(3);
surf(X, T, Difft1);
xlabel('x');
ylabel('t');
zlabel('$|\Phi_D(x,t)*N(x)-u_1(x,t)|$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca, 'hot');
title('$|\Phi_D(x,t)*N(x)-u_1(x,t)|$', 'Interpreter', 'latex');

subplot(2,4,7);
view(3);
surf(X, T, Difft2);
xlabel('x');
ylabel('t');
zlabel('$|\Phi_D(x,t)*N(x)-u_2(x,t)|$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap(gca, 'hot');
title('$|\Phi_D(x,t)*N(x)-u_2(x,t)|$', 'Interpreter', 'latex');
</source>

= Referencias =

* Sandro Salsa, Gianmaria Verzini. Partial Differential Equations in Action From Modelling to Theory ([https://link.springer.com/book/10.1007/978-3-031-21853-8 [1<nowiki>]</nowiki>]).

[[Categoría:EDP]]
[[Categoría:EDP24/25]]

Ecuación del calor (ADMR)

2025-03-18T22:56:44Z

Daniel Rodríguez Calderón: /* Códigos */

{{ TrabajoED | Ecuación del calor (Grupo ADMR). | [[:Categoría:EDP|EDP]]|[[:Categoría:EDP24/25|2024-25]] |
*Ángel De Lucas Miranda,
*Daniel Rodríguez Calderón,
*Marcos Cabellos Hernández,
*Rafael Pascual Ortega.}}

=Introducción=

El desarrollo de la ecuación del calor se atribuye principalmente al matemático y físico francés Joseph Fourier, quien en 1822 publicó su obra ''Théorie analytique de la chaleur''. En este trabajo formuló una ecuación diferencial en derivadas parciales que describe la propagación del calor por difusión en un medio continuo. En esta ocasión, profundicemos en dos aspectos de este modelo.

Desde un punto de vista físico, estas ecuaciones permiten la modelización de fenómenos como el calentamiento de materiales o la difusión de solutos en fluidos estacionarios. Sin embargo, la influencia de la física es más profunda: de hecho, el coeficiente de difusión lleva su nombre por la interpretación física de su acción en el sistema. ¿Cuál es este efecto, y por qué se lo denota de esta forma?

Planteemos también el paso entre problemas de dominio acotado y los de no acotado. ¿Podemos aproximar las soluciones del sistema no acotado mediante las del acotado? Y de ser así, ¿qué condiciones de frontera en una dimensión funcionan mejor?

=Solución acotada vs Solución no acotada=
Demos por conocidos el uso de separación de variables para resolver problemas de autofunciones en caso acotado y el manejo de la solución fundamental junto con la convolución para el sistema no acotado. Para comparar todas las soluciones correspondientes, consideremos sistema <math> 1-</math>dimensional, con una gaussiana <math> e^{-x^2} </math> como condición inicial. Justificamos físicamente como el uso de un láser para calentar un punto de una tira fina de metal. Considerando que los fotones puedan desviarse como una distribución normal por el teorema central del límite, es un montaje experimental factible y regular. ¿Y qué condiciones de frontera?

Si consideramos de Dirichlet, a una tira de longitud <math> 2a </math> le colocaremos en los extremos un material que se mantenga siempre a la misma temperatura, pongamos que a una temperatura nula. El sistema resulta

<math>\quad
\begin{align}
\begin{cases} u_{t} - u_{xx} = 0 \quad &\forall t > 0, x\in [-a,a],\\
u(-a,t) = u(a,t) = 0 \quad &\forall t > 0, \\
u(x,0) = e^{-x^2} \quad &\forall x\in [-a,a].
\end{cases}
\end{align}
</math>

La solución al resolver el problema de autofunciones viene dada por

<math>\quad
u(x,t) = \sum_{n=1}^\infty A_{n} \cos\left(\frac{\pi}{a}(n - \frac{1}{2}) x\right) e^{-\dfrac{\pi^2}{a^2}\left(n - \frac{1}{2}\right)^2 t},
\quad \text{ donde } \quad
A_n = \frac{\int_{-a}^{a} u(x,0) \cos\left(\frac{\pi}{a}(n - \frac{1}{2}) x\right)dx}{\int_{-a}^{a} \cos^2\left(\frac{\pi}{a}(n - \frac{1}{2}) x\right)dx}.
</math>

Consideremos ahora condiciones de frontera Neumann, con un flujo nulo que físicamente viene de poner extremos aislantes. Para cierto <math> b > 0,</math> y no necesariamente igual que con <math> a </math>. Consideraremos la ecuación del calor

<math>\quad
\begin{align}
\begin{cases} u_{t} - u_{xx} = 0 \quad &\forall t > 0 , x\in [-b,b],\\
u_x(-b,t) = u_x(b,t) = 0 \quad &\forall t > 0, \\
u(x,0) = e^{-x^2} \quad &\forall x\in[-b,b].
\end{cases}
\end{align}
</math>

De nuevo podemos resolver el sistema, por ser homogéneo, mediante separación de variables, obteniendo

<math>\quad
u(x,t) = \sum_{m=0}^\infty B_{m} cos\left(\frac{m\pi}{b} x\right) e^{-\dfrac{m^2\pi^2}{b^2} t},\quad \text{ donde } \quad
B_m = \frac{\int_{-b}^{b} u(x,0) \cos\left(\frac{m\pi}{b} x\right)dx}{\int_{-b}^{b} \cos^2\left(\frac{m\pi}{b}x\right)dx}.
</math>

Un último detalle antes de proceder a los errores es la obtención simbólica de la solución no acotada.

<math>\quad
\begin{align}
\begin{cases} u_{t} - u_{xx} = 0 \quad &\forall t>0, x\in \mathbb{R},\\
u(x,0) = e^{-x^2} \quad &\forall x\in \mathbb{R}
\end{cases}
\end{align}
\quad
\Longrightarrow
\quad
u(x,t) = \frac{1}{\sqrt{4\pi t}}e^{\frac{-x^2}{4t}} * e^{-x^2} = \frac{e^{\frac{-x^2}{4t+1}}}{\sqrt{4t+1}}.
</math>

Queremos comprobar cuán bien aproximan las soluciones del caso acotado Dirichlet y Neumann a la obtenida por la solución fundamental. Para ello, escojamos un tiempo máximo <math>t_{max} </math> hasta el cuál tomaremos el error en norma <math> L^{\infty} </math>. Queremos estudiar:

* El efecto de aumentar el rango de definición con <math> a,b </math>.
* El efecto de aumentar el número de términos en las series de los sistemas acotados.

=== Análisis de los errores ===
[[Archivo:GrupoADMR Calor Errores1.jpg|550px|thumb|right| Representación gráfica de las soluciones del problema acotado (<math>a = b = 20</math>) y no acotado truncando las series de Fourier en su término <math>100</math>, para un tiempo máximo de <math>1</math>; y de sus errores. Nótese el orden de magnitud de los errores -pequeño-, los picos generados por los errores numéricos dados por la fórmula del trapecio y la suavidad en la frontera dada por las condiciones en esta.]]

Consideramos los parámetros, donde <math> n </math> es el término de la serie de Fourier donde truncamos:

*<math> a = b = 20, </math>
*<math> t_{max} = 1, </math>
*<math> n = 100. </math>

Con <math> t_{max} </math> fijo podemos observar que al aumentar <math>a, b</math> y <math>n</math> el error entre las soluciones acotadas y no acotada disminuye. Nótese que se estanca a partir de cierto <math>n</math>. De hecho, el error llega a aumentar a partir de ciertos <math> a,b </math> si truncamos demasiado pronto o demasiado tarde por la aproximación del trapecio para particiones del espacio demasiado bastas -no observado por pericia en código final-.

{|style="margin: 0 auto;"
| [[Archivo:GrupoADMR Calor Errores MaxD.jpeg|250px|thumb|upright|]]
| [[Archivo:GrupoADMR Calor Errores MaxNe.jpeg|250px|thumb|upright|]]
|}

Al ampliar el tiempo con la <math>n</math> y fijada la longitud, al hacer más grande <math>t_{max}</math> los errores máximos no aumentan. En cambio, si tomamos el máximo error para diferentes tiempos vemos con suficiente precisión que no hay diferencia pero que para truncamientos precoces Neumann tiene ligeramente menor error.

{|style="margin: 0 auto;"
| [[Archivo:GrupoADMR Calor Errores DNo.jpeg|250px|thumb|upright|]]
| [[Archivo:GrupoADMR Calor Errores NoNe.jpeg|250px|thumb|upright|]]
| [[Archivo:GrupoADMR Calor Errores DNe.jpeg|250px|thumb|upright|]]
|}

=Coeficiente de difusión=

[[Archivo:Fundamental.gif|300px|thumb|right|Solución fundamental de la ecuación del calor unidimensional para distintos valores del coeficiente de difusión. La representación gráfica parte del tiempo inicial <math/> 0.01, <math/> esperando una delta de Dirac en el instante inicial. Nótese que en este tiempo, a mayor coeficiente de difusión más se ha difundido el calor en el sistema, lo que se nota especialmente en el valor máximo observado y su disminución cuanto más crece el coeficiente.]]

Vamos a interpretar como influye exactamente el coeficiente de difusión en la solución de la ecuación del calor. Imaginemos por un segundo que tenemos varias barras (cada una de un material diferente) y queremos estudiar como se difunde el calor a través de cada una de ellas. Usando Matlab podemos calcular y representar fácilmente la solución del problema para diferentes valores del coeficiente de difusión, que denotaremos por <math/> D <math/>.

## Interpretemos el valor de <math/> D <math/> aprovechando el montaje físico anterior. Numéricamente vemos que cuanto mayor es este coeficiente más rápido se difunde el calor, evidenciado en la caída de la solución en tiempos cada vez más cortos.

[[Archivo:Fundamental.gif|400px|thumb|right|Solución fundamental de la ecuación del calor unidimensional para distintos valores del coeficiente de difusión. La representación gráfica parte del tiempo inicial <math/> 0.01, <math/> esperando una delta de Dirac en el instante inicial. Nótese que en este tiempo, a mayor coeficiente de difusión más se ha difundido el calor en el sistema, lo que se nota especialmente en el valor máximo observado y su disminución cuanto más crece el coeficiente.]]

Como se puede observar en la animación, fijado un tiempo, a medida que crece el valor del coeficiente de difusión el pico de la solución fundamental se va haciendo más cada vez más pequeño, es decir, a mayor coeficiente de difusión, menor es el tiempo que tarda en expandirse el calor. De hecho, podemos observar este hecho directamente de la expresión de la solución fundamental del calor. La solución de la ecuación del calor unidimensional en una barra de longitud $ L = 1 $, con coeficiente de difusión $ D $ y condición inicial dada por una función gaussiana $ f(x) = e^{-x^2} $, se expresa como una serie de Fourier:

## En el caso unidimensional acotado con solución

<math>
u(x,t) = \sum_{n=1}^{\infty} A_n e^{-D (n\pi)^2 t} \sin(n\pi x),
</math>

[Solución que depende de D]- [Viendo la solución del apartado 3 me estoy dando cuenta que también depende de la longitud de la barra, si estuviésemos trabajando con el problema acotado podríamos ver como afectan estos dos parámetros al tiempo de difusión]

Tomando el límite de D tendiendo a <math/> \infty <math/> vemos que u(x,t) tiende a 0, fijado un t determinado. Por tanto, podemos concluir que la solución del calor alcanzará el estado estacionario antes si su coeficiente del calor es mayor.

## Podemos ver que un mayor coeficiente <math/> D <math/> acelera el decaimiento de la exponencial por lo que tendrá un efecto similar justificando su nombre.

= Códigos =

<source lang="matlab">
%%% Código generado junto con Chat GPT para la representación de la solución fundamental de la ecuación del calor unidimensional para distintos valores del coeficiente de difusión

clc
clear all
close all

% Vector de valores de D (coeficiente de difusión)
DD = [1,2,3,4,5,6,7,8,9,10];

% Nombre del archivo GIF de salida
output_gif = 'Fundamental.gif';

% Crear la figura
figura = figure(1);
grid on;
view(3);

% Bucle para iterar sobre cada valor de D
for i = 1:length(DD)
D = DD(i);

% Función fundamental de la ecuación del calor
phi_D = @(x, t) exp(-x.^2./(4*D*t))./sqrt(4*D*pi*t);

% Definición de la malla
xx = -1:0.01:1; % Rango de x
tt = 10^(-2):0.01:1; % Rango de t
[X, T] = meshgrid(xx, tt); % Malla de coordenadas
PHI_D = phi_D(X, T); % Evaluación de la función

% Graficar la superficie
clf;
surf(X, T, PHI_D);
shading flat
shading interp
colormap('jet')
xlabel('x');
ylabel('t');
zlabel('$\Phi_D(x,t)$', 'Interpreter', 'latex');
title(['Solución fundamental de la ecuación del calor, cuando D=', num2str(D)])
axis([-1, 1, 0, 1, 0, 3])
drawnow;

% Captura de la imagen para el GIF
frame = getframe(figura);
img = frame2im(frame);
[imind, cm] = rgb2ind(img, 256);

% Escribir la imagen en el archivo GIF
if i == 1
imwrite(imind, cm, output_gif, 'gif', 'LoopCount', Inf, 'DelayTime', 0.7);
else
imwrite(imind, cm, output_gif, 'gif', 'WriteMode', 'append', 'DelayTime', 0.7);
end

end
</source>

[[Archivo:Convoluciones_teorica.gif|500px|thumb|right|Soluciones de la ecuación del calor unidimensional para distintos valores del coeficiente de difusión.]]

<source lang="matlab">
%%% Código generado junto con Chat GPT para la representación de la soluciones de la ecuación del calor unidimensional para distintos valores del coeficiente de difusión

clc
clear all
close all

% Vector de valores de D (coeficiente de difusión)
DD = [1,10,20,30,40,50,60,70,80,90,100];

% Función gaussiana
gaus = @(x) (exp(-x.^2));

% Nombre del archivo GIF de salida
output_gif = 'Convoluciones_teorica.gif';

% Crear la figura
figura = figure(1);
grid on;
view(3);

% Bucle para iterar sobre cada valor de D
for h = 1:length(DD)
D = DD(h);

% Funciones fundamentales y convolución
phi_D = @(x, t) exp(-x.^2./(4*D*t))./sqrt(4*D*pi*t);
convo = @(x,t) (exp(-x.^2./(4*D.*t+1))/sqrt(4*D.*t+1));

% Definición de la malla
xx = -10:0.1:10; % Rango de x
tt = 10^(-2):0.01:1; % Rango de t

% Inicialización de la matriz de convolución
conv = zeros(length(xx), length(tt));

% Cálculo de la convolución
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
val_conv(i,j) = convo(auxx, auxt);
end
end

% Graficar la superficie
clf;
hold on;
view(3);
surf(xx, tt, val_conv.');
xlabel('x');
ylabel('t');
zlabel('$\Phi_D(x,t)*N(x)$', 'Interpreter', 'latex');
shading flat;
shading interp;
colormap('jet');
title(['Soluciones en dominio no acotado, cuando D=', num2str(D)]);
axis([-10, 10, 0, 1, 0, 1]);
hold off;
drawnow;

% Captura de la imagen para el GIF
frame = getframe(figura);
img = frame2im(frame);
[imind, cm] = rgb2ind(img, 256);

% Escribir la imagen en el archivo GIF
if h == 1
imwrite(imind, cm, output_gif, 'gif', 'LoopCount', Inf, 'DelayTime', 0.7);
else
imwrite(imind, cm, output_gif, 'gif', 'WriteMode', 'append', 'DelayTime', 0.7);
end

end
</source>

<source lang="matlab">
%%% Código generado junto con Chat GPT para la representación de los errores entre las soluciones

function [Diff12,Difft1,Difft2] = errores2(n_max, a, b,t_max)

% Función como condición inicial: una función Gaussiana
gaus = @(x) (exp(-x.^2));
u0 = @(x) exp(-x.^2);
D = 1; % Coeficiente de difusión

% Funciones fundamentales y convolución para el modelo de difusión
phi_D = @(x, t) exp(-x.^2./(4*D*t))./sqrt(4*D*pi*t);
convo = @(x,t) (exp(-x.^2./(4*D.*t+1))/sqrt(4*D.*t+1));

% Definición de la malla (espacio y tiempo)
xx = -a:0.01:a;
tt = 10^(-2):0.01:t_max;
[X, T] = meshgrid(xx, tt);

% Inicialización de la matriz de convolución
conv = zeros(length(xx), length(tt));

% Cálculo de la convolución
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
val_conv(j,i) = convo(auxx, auxt);
end
end

% Inicialización de la matriz de resultados para la solución Dirichlet
res_dir = zeros(length(xx), length(tt));

% Cálculo de la solución directa
for n=1:n_max
A(n) = trapz(xx, u0(xx) .* cos((pi/a)*(n - 1/2) * xx)) /...
trapz(xx, cos((pi/a)*(n - 1/2) * xx).^2);
end

% Evaluación de la solución directa en el dominio espacio-temporal
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
u = 0;

% Sumar los términos de la serie de Fourier
for n=1:n_max
u = u + A(n) .* cos((pi/a)*(n - 1/2) * auxx) .*...
exp(-((pi^2 / a^2) * (n - 1/2)^2) .* auxt);
end
res_dir(i,j) = u;
end
end

% Inicialización de la matriz de resultados para la solución en frontera Neumann
res_neu = zeros(length(xx), length(tt));

% Cálculo de la solución en frontera Neumann
for n=1:n_max
B(n) = trapz(xx, u0(xx) .* cos((n*pi/b) * xx)) /...
trapz(xx, cos((n*pi/b) * xx).^2);
end

% Evaluación de la solución en frontera Neumann en el dominio espacio-temporal
for j = 1:length(tt)
auxt = tt(j);
for i = 1:length(xx)
auxx = xx(i);
u = trapz(xx, u0(xx))./(2.*b);

% Sumar los términos de la serie de Fourier
for n=1:n_max
u = u + B(n) .* cos((n*pi/b) * auxx) .* exp(-((n^2 * pi^2 / b^2)) .* auxt);
end
res_neu(i,j) = u;
end
end

% Transponer los resultados para compararlos con la convolución
U1 = res_dir.';
U2 = res_neu.';

% Calcular las diferencias entre las soluciones y la convolución
Diff12 = abs(U1 - U2);
Difft1 = abs(val_conv - U1);
Difft2 = abs(val_conv - U2);

end
</source>

= Referencias =

* Sandro Salsa, Gianmaria Verzini. Partial Differential Equations in Action From Modelling to Theory ([https://link.springer.com/book/10.1007/978-3-031-21853-8 [1<nowiki>]</nowiki>]).

[[Categoría:EDP]]
[[Categoría:EDP24/25]]

Ecuación del calor (ADMR)

2025-03-18T22:39:43Z

Daniel Rodríguez Calderón: /* Coeficiente de difusión */

Ecuación del calor (ADMR)

2025-03-18T22:32:58Z

Daniel Rodríguez Calderón: /* Análisis de los errores */

Ecuación del calor (ADMR)

2025-03-18T22:32:32Z

Daniel Rodríguez Calderón: /* Solución acotada vs Solución no acotada */

Ecuación del calor (ADMR)

2025-03-18T22:23:29Z

Daniel Rodríguez Calderón: /* Introducción */

Ecuación del calor (ADMR)

2025-03-18T22:06:23Z

Daniel Rodríguez Calderón: /* Solución acotada vs Solución no acotada */

Ecuación del calor (ADMR)

2025-03-18T22:04:50Z

Daniel Rodríguez Calderón: /* Coeficiente de difusión */

Ecuación del calor (ADMR)

2025-03-18T22:04:13Z

Daniel Rodríguez Calderón: /* Solución acotada vs Solución no acotada */

Ecuación del calor (ADMR)

2025-03-16T16:49:26Z

Daniel Rodríguez Calderón: /* Códigos */

Ecuación del calor (ADMR)

2025-03-16T16:49:02Z

Daniel Rodríguez Calderón: /* Códigos */

Ecuación del calor (ADMR)

2025-03-16T16:48:45Z

Daniel Rodríguez Calderón: /* Códigos */

Ecuación del calor (ADMR)

2025-03-16T16:48:25Z

Daniel Rodríguez Calderón: /* Códigos */

Ecuación del calor (ADMR)

2025-03-16T16:47:17Z

Daniel Rodríguez Calderón: /* Códigos */

Ecuación del calor (ADMR)

2025-03-16T16:23:18Z

Daniel Rodríguez Calderón: /* Solución acotada vs Solución no acotada */

Ecuación del calor (ADMR)

2025-03-16T16:22:49Z

Daniel Rodríguez Calderón: /* Solución acotada vs Solución no acotada */

Ecuación del calor (ADMR)

2025-03-16T16:22:15Z

Daniel Rodríguez Calderón: /* Códigos */

Ecuación del calor (ADMR)

2025-03-16T16:21:18Z

Daniel Rodríguez Calderón: /* Códigos */

Ecuación del calor (ADMR)

2025-03-16T16:20:47Z

Daniel Rodríguez Calderón: /* Códigos */

Ecuación del calor (ADMR)

2025-03-16T16:16:38Z

Daniel Rodríguez Calderón: /* Códigos */

Archivo:Convoluciones teorica.gif

2025-03-16T16:15:51Z

Daniel Rodríguez Calderón: Soluciones de la ecuación del calor para diferentes coeficientes de disfusión.

Soluciones de la ecuación del calor para diferentes coeficientes de disfusión.

Ecuación del calor (ADMR)

2025-03-16T16:14:35Z

Daniel Rodríguez Calderón: /* Códigos */

Ecuación del calor (ADMR)

2025-03-16T16:14:18Z

Daniel Rodríguez Calderón: /* Códigos */

Ecuación del calor (ADMR)

2025-03-16T16:13:30Z

Daniel Rodríguez Calderón: /* Códigos */

Ecuación del calor (ADMR)

2025-03-16T15:26:12Z

Daniel Rodríguez Calderón: /* Códigos */

Ecuación del calor (ADMR)

2025-03-16T15:23:45Z

Daniel Rodríguez Calderón: /* Códigos */

Archivo:Fundamental.gif

2025-03-16T15:22:33Z

Daniel Rodríguez Calderón: Solución fundamental para la ecuación del calor con distintos valores de coeficiente de difusión.

Solución fundamental para la ecuación del calor con distintos valores de coeficiente de difusión.