Diferencia entre revisiones de «Ecuación del calor (Grupo ACIRV)»

De MateWiki
Saltar a: navegación, buscar
(→‎Modelización de la ecuación del calor con una dimensión)
(→‎Comparación de resultados y análisis de errores)
 
(No se muestran 45 ediciones intermedias de 4 usuarios)
Línea 3: Línea 3:
  
 
=Introducción=
 
=Introducción=
 +
 +
En esta página ilustraremos el uso del método de diferencias finitas para resolver la ecuación del calor en una dimensión.
 +
 +
Dicho método sirve para resolver de manera aproximada ecuaciones diferenciales parciales y ordinarias de manera numérica. Funciona sustituyendo las derivadas por cocientes de diferencias que vienen dados por valores de la función en puntos discretos de una malla.
 +
  
  
 
=Modelización de la ecuación del calor con una dimensión=
 
=Modelización de la ecuación del calor con una dimensión=
  
−
En primer lugar vamos a resolver el siguiente problema de manera analítica, para después comparar su solución con la obtenida numéricamente.
+
La ecuación a resolver es
  
−
Se considera una varilla metálica que ocupa el intervalo \([0,1]\) y que se encuentra aislada por su superficie lateral, de manera que la conducción de calor solo se produce en la dirección longitudinal. La temperatura inicial de la varilla es \(10^{\circ}C\). En el extremo derecho se consigue mantener la temperatura a \(10^{\circ}C\) mientras que en el izquierdo la temperatura es siempre de \(1^{\circ}C\).
+
<center><math>
 +
u_t = \alpha u_{xx}, \quad x \in [0,1], \ t > 0,\ \alpha = 1.
 +
</math></center>
  
 +
donde \(u(x,t)\) representa la temperatura en el punto \(x\) y en el instante \(t\), y \(\alpha\) es una constante que viene dada por la conductividad térmica del material.
 +
 +
Se considera una varilla metálica en el intervalo \([0,1]\) y aislada en su superficie lateral tal que la conducción de calor solo se produce en la dirección longitudinal. La temperatura inicial de la varilla es \(10^{\circ}C\). En los extremos derecho e izquierdo la temperatura se mantiene a \(10^{\circ}C\)  y \(1^{\circ}C\) respectivamente.
  
 
El sistema que modeliza el comportamiento de la temperatura, representada por la función \(u(x,t)\), es el siguiente:
 
El sistema que modeliza el comportamiento de la temperatura, representada por la función \(u(x,t)\), es el siguiente:
Línea 17: Línea 27:
 
\left \{
 
\left \{
 
\begin{array}{llll}
 
\begin{array}{llll}
−
u_t - u_{xx} & = 0 & 0<x<1, & t>0 \\
+
u_t - u_{xx} & = 0 & x \in [0,1], & t>0 \\
 
u(0,t)& = 1 & t>0\\
 
u(0,t)& = 1 & t>0\\
 
u(1,t) & = 10 & t>0\\
 
u(1,t) & = 10 & t>0\\
−
u(x,0) &= 10 & 0<x<1
+
u(x,0) &= 10 & x \in [0,1]
 
\end{array}
 
\end{array}
 
\right .
 
\right .
 +
\tag{1}
 
</math></center>
 
</math></center>
  
−
Puesto que no es un sistema homogéneo, debemos encontrar primero la solución estacionaria. Es decir, suponemos que \(t \rightarrow \infty\) (la solución ya no varía en el tiempo).
+
=Resolución analítica del problema=
  
−
La solución obtenida es: \(v(x) = 9x +1\)
+
Primeramente, resolveremos el problema analíticamente para después comparar su solución con la obtenida numéricamente.
  
−
Una vez calculada la solución estacionaria, homogeneizamos el sistema definiendo la ecuación \(w(x,t) = u(x,t) - v(x)$.
+
Puesto que el sistema obtenido \((1)\) no es homogéneo, debemos encontrar primero la solución estacionaria. Es decir, suponemos que \(t \rightarrow \infty\) (la solución ya no varía en el tiempo).
 +
 
 +
La solución estacionaria obtenida es \(v(x) = 9x +1\) . Al representarla gráficamente en 3D obtenemos:
 +
 
 +
[[Archivo:Sol estacionaria.png|600px|thumb|right|Solución estacionaria]]
 +
 
 +
{{matlab|codigo=%Creamos una matriz que representa una malla de puntos (tiempo,
 +
% espacio) y en cada una de las columnas introducimos el valor de 9x+1 para todos los tiempos.
 +
 
 +
% Matriz
 +
evaluaciones = zeros(100,100);
 +
 
 +
% Límite temporal y vectores
 +
T = 1;
 +
t = linspace(0,T,100);
 +
x = linspace(0,1,100);
 +
 
 +
% Rellenamos la matriz
 +
for j = 1:100
 +
    evaluaciones(:,j) = (9*x(j) + 1) * ones(100,1);
 +
end
 +
 
 +
% Representación gráfica
 +
surf(t,x,evaluaciones')
 +
title('Solución estacionaria: v(x) = 9x + 1')
 +
xlabel('Tiempo')
 +
ylabel('Espacio')
 +
zlabel('Temperatura (°C)')
 +
 
 +
}}
 +
 
 +
Una vez calculada la solución estacionaria, homogeneizamos el sistema definiendo la ecuación \(w(x,t) = u(x,t) - v(x)\).
  
 
<center><math>
 
<center><math>
 
\left \{
 
\left \{
 
\begin{array}{llll}
 
\begin{array}{llll}
−
w_t - w_{xx} & = 0 & 0<x<1, & t>0 \\
+
w_t - w_{xx} & = 0 & x \in [0,1], & t>0 \\
 
w(0,t)& = 0 & t>0\\
 
w(0,t)& = 0 & t>0\\
 
w(1,t) & = 0 & t>0\\
 
w(1,t) & = 0 & t>0\\
−
w(x,0) &= 9(1-x) & 0<x<1
+
w(x,0) &= 9(1-x) & x \in [0,1]
 
\end{array}
 
\end{array}
 
\right .
 
\right .
Línea 44: Línea 86:
 
Resolviendo mediante el método de separación de variables y por superposición de soluciones, obtenemos la siguiente solución:
 
Resolviendo mediante el método de separación de variables y por superposición de soluciones, obtenemos la siguiente solución:
  
 +
<center><math>
 +
w(x,t) = \sum^{\infty}_{k=1} \frac{18}{k\pi}sen(k\pi x)e^{-k^2 \pi^2 t} .
 +
</math></center>
 +
 +
Véase la evolución de la temperatura en función del tiempo en cada punto de la barra. Representamos la temperatura respecto de la posición, variando el tiempo hasta \(t = 0.5 s\). Podemos representarla para diferentes términos, veamos que utilizando 100 se aproxima suficiente a tomar infinitos términos.
 +
 +
 +
[[Archivo:Heat equation analitico.gif|600px|thumb|right|Solución analítica]]
 +
[[Archivo:Heat equation analitico 3.gif|400px|thumb|right|Solución analítica 2 términos]]
 +
[[Archivo:Heat equation analitico 2.gif|400px|thumb|right|Solución analítica 10 términos]]
 +
 +
{{matlab|codigo=
 +
% Resolución de la ecuación del calor 1D por diferencias finitas y comparación con solución analítica
 +
clear all; close all;
 +
% Parámetros
 +
L = 1;        % Longitud de la varilla
 +
Nx = 20;      % Número de puntos en el espacio
 +
dx = L / (Nx-1);
 +
alpha = 1; % Difusividad térmica
 +
dt = 0.0005;  % Paso de tiempo
 +
Nt = 1000;    % Número de pasos de tiempo
 +
 +
% Discretización espacial
 +
x = linspace(0, L, Nx);
 +
 +
% Condiciones iniciales
 +
u = 10 * ones(Nx, 1);
 +
 +
% Condiciones de frontera
 +
u(1) = 1;  % Extremo izquierdo
 +
u(end) = 10; % Extremo derecho
 +
 +
% Definimos la solucón estacionaria
 +
 +
estacionaria = @(x)9*x +1;
 +
 +
% Coeficiente del método
 +
r = alpha * dt / dx^2;
 +
if r > 0.5
 +
    error('El método no es estable, reduce dt o aumenta dx.');
 +
end
 +
 +
% Preparar la grabación de un GIF
 +
filename = 'heat_equation_analitico.gif';
 +
frame_delay = 0.1;
 +
 +
% Evolución temporal
 +
for n = 1:Nt
 +
    % Solución analítica
 +
    t = n * dt;
 +
    w = zeros(size(x));
 +
    for k = 1:100 % Consideramos los primeros 100 términos
 +
%        w = w + (18 / (k * pi)) * sin(k * pi * x) * exp(-k^2 * pi^2 * t);
 +
        w = w + (-(18*(sin(pi*k) - pi*k))/(k^2*pi^2))*sin(k * pi * x) * exp(-k^2 * pi^2 * t);
 +
 +
    end
 +
   
 +
    % Guardar frames en el GIF
 +
    if mod(n, 10) == 0
 +
        plot(x, estacionaria(x), 'g-', 'LineWidth', 2);
 +
        hold on;
 +
        plot(x, w + 9*x+ 1, 'r--', 'LineWidth', 2); % +1 para respetar la condición de frontera izquierda
 +
        hold off;
 +
        xlabel('Posición'); ylabel('Temperatura');
 +
        title(sprintf('Tiempo: %.3f s', n*dt));
 +
        legend('Solución estacionaria','Solución Analítica');
 +
        drawnow;
 +
        frame = getframe(gcf);
 +
        img = frame2im(frame);
 +
        [A, map] = rgb2ind(img, 256);
 +
        if n == 10
 +
            imwrite(A, map, filename, 'gif', 'LoopCount', inf, 'DelayTime', frame_delay);
 +
        else
 +
            imwrite(A, map, filename, 'gif', 'WriteMode', 'append', 'DelayTime', frame_delay);
 +
        end
 +
    end
 +
end
 +
 +
% Gráfica final
 +
figure;
 +
plot(x, estacionaria(x), 'g-', 'LineWidth', 2);
 +
hold on;
 +
plot(x, w + 9*x + 1, 'r--', 'LineWidth', 2);
 +
hold off;
 +
xlabel('Posición'); ylabel('Temperatura');
 +
title('Comparación final de temperatura en la varilla');
 +
legend('Solución estacionaria', 'Solución Analítica');
 +
}}
 +
 +
Vemos que según aumenta el tiempo, se ajusta a la solución estacionaria.
 +
 +
=Resolución numérica de la ecuación del calor=
 +
 +
Falta obtener la resolución numérica del sistema planteado anteriormente \((1)\).
 +
 +
Así, discretizaremos el espacio en una malla de puntos espaciales \( x_i \) y tiempos \( t_n \). Estos deben cumplir \(x_0 = 0, \ x_N = 1\) de acuerdo con los intervalos de definición de la EDP. Definimos \(\Delta t = t_{n+1} - t_n\) y \(\Delta x = x_{i+1} - x_i\).
 +
 +
Pasamos a discretizar la ecuación según el método de diferencias finitas explícito. De \(u_t = \alpha u_{xx}\), pasamos a
 +
 +
<center><math>
 +
\frac{u_{i}^{n+1} - u_{i}^{n}}{\Delta t} = \alpha \frac{u_{i+1}^{n} - 2u_{i}^{n} + u_{i-1}^{n}}{\Delta x^2}.
 +
</math></center>
 +
 +
Despejando \( u_{i}^{n+1} \), obtenemos la ecuación de actualización
 +
 +
<center><math>
 +
u_{i}^{n+1} = u_{i}^{n} + r \left( u_{i+1}^{n} - 2 u_{i}^{n} + u_{i-1}^{n}\right) .
 +
</math></center>
 +
 +
Definimos el número de Fourier
 +
\(r = \frac{\alpha \Delta t}{\Delta x^2}
 +
\), que mide la relación entre la difusión térmica, el paso de tiempo y el desplazamiento. En nuestro caso tenemos \(\alpha = 1\).
 +
 +
 +
Por último, inicializamos el problema para las condiciones de contorno. Definimos \( u_0^n = 1 , u_N^n = 10 \). Y las condiciones iniciales \(u_i^0 = 10, \ \forall i = 1,...,N.\)
 +
 +
Vamos a representar la solución numérica para diferentes mallados
 +
 +
 +
[[Archivo:Heat equation num 4.gif|600px|thumb|right|Solución numérica para N = 3 (4 puntos y \(\Delta x = 0.33\)), \(\Delta t = 0.0005\) ]]
 +
[[Archivo:Heat equation num 20.gif|600px|thumb|right|Solución numérica para N = 20 (21 puntos y \(\Delta x = 0.05\)), \(\Delta t = 0.0005\)]]
 +
 +
{{matlab|codigo=
 +
% Resolución de la ecuación del calor 1D por diferencias finitas
 +
 +
% Parámetros
 +
L = 1;        % Longitud de la varilla
 +
Nx = 21;      % Número de puntos en el espacio
 +
dx = L / (Nx-1);
 +
D = 1; % Difusividad térmica
 +
dt = 0.0005;  % Paso de tiempo
 +
Nt = 1000;    % Número de pasos de tiempo
 +
 +
% Discretización espacial
 +
x = linspace(0, L, Nx);
 +
 +
% Condiciones iniciales
 +
u = 10 * ones(Nx, 1);
 +
 +
% Condiciones de frontera
 +
u(1) = 1;  % Extremo izquierdo
 +
u(end) = 10; % Extremo derecho
 +
 +
%Solución estacionaria
 +
estacionaria = @(x)9*x +1;
 +
 +
 +
% Coeficiente del método
 +
r = D * dt / dx^2;
 +
if r > 0.5
 +
    error('El método no es estable, reduce dt o aumenta dx.');
 +
end
 +
 +
% Preparar la grabación de un GIF
 +
filename = 'heat_equation_num_20.gif';
 +
frame_delay = 0.1;
 +
 +
% Evolución temporal
 +
for n = 1:Nt
 +
    u_new = u; % Copia de la solución actual
 +
    for i = 2:Nx-1
 +
        u_new(i) = u(i) + r * (u(i+1) - 2*u(i) + u(i-1));
 +
    end
 +
    u = u_new;
 +
   
 +
    % Guardar frames en el GIF
 +
    if mod(n, 10) == 0
 +
        plot(x, u, 'b-', 'LineWidth', 2);
 +
        hold on;
 +
        plot(xx, estacionaria(xx), 'g-', 'LineWidth', 2); % Solución estacionaria
 +
        hold off;        xlabel('Posición'); ylabel('Temperatura');
 +
        title(sprintf('Tiempo: %.3f s', n*dt));
 +
        legend('Diferencias Finitas', 'Solución estacionaria');
 +
        drawnow;
 +
        frame = getframe(gcf);
 +
        img = frame2im(frame);
 +
        [A, map] = rgb2ind(img, 256);
 +
        if n == 10
 +
            imwrite(A, map, filename, 'gif', 'LoopCount', inf, 'DelayTime', frame_delay);
 +
        else
 +
            imwrite(A, map, filename, 'gif', 'WriteMode', 'append', 'DelayTime', frame_delay);
 +
        end
 +
    end
 +
end
 +
 +
% Gráfica final
 +
figure;
 +
plot(x, u, 'b-', 'LineWidth', 2);
 +
xlabel('Posición'); ylabel('Temperatura');
 +
 +
title('Distribución final de temperatura en la varilla');
 +
 +
}}
 +
 +
=Análisis de estabilidad=
 +
 +
Es importante asegurarse de que el método es estable. Aplicamos el análisis de von Neumann, en el que asumimos soluciones de la forma:
 +
 +
<center><math>
 +
u_{i}^{n} = G^n e^{ikx_i} ,
 +
</math></center>
 +
 +
donde \( G \) es el factor de amplificación. Reemplazando esta forma en la ecuación de actualización y tras un análisis algebraico, se obtiene:
 +
 +
<center><math>
 +
G = 1 - 4r \sin^2 \left( \frac{k \Delta x}{2} \right) .
 +
</math></center>
 +
 +
Para que el método sea estable, es necesario que \( |G| \leq 1 \). Como el término \( \sin^2 (\cdot) \) está acotado entre 0 y 1, se deduce la condición de estabilidad:
 +
 +
<center><math>
 +
1 - 4r \leq 1 \Rightarrow r \leq \frac{1}{2} .
 +
</math></center>
 +
 +
Si \( r > \frac{1}{2} \), entonces \( |G| > 1 \), haciendo que la solución sea inestable.
 +
 +
=Comparación de resultados y análisis de errores=
 +
 +
Seguidamente, compararemos la solución analítica con la numérica para verificar la precisión de esta última en función de la malla empleada.
 +
 +
La siguiente representación recoge cómo varía la temperatura según qué solución en función del tiempo en cada punto.
 +
 +
[[Archivo:Heat equation comparacion 2.gif|600px|thumb|right|Comparación solución analítica y numérica para 4 puntos]]
 +
 +
{{matlab|codigo=
 +
% Resolución de la ecuación del calor 1D por diferencias finitas y comparación con solución analítica
 +
 +
% Parámetros
 +
L = 1;        % Longitud de la varilla
 +
Nx = 4;      % Número de puntos en el espacio
 +
dx = L / (Nx-1);
 +
alpha = 1; % Difusividad térmica
 +
dt = 0.0005;  % Paso de tiempo
 +
Nt = 1000;    % Número de pasos de tiempo
 +
 +
% Discretización espacial
 +
x = linspace(0, L, Nx); %Discretización para hacer los cálculos con la numérica
 +
xx = linspace(0,L,50); %Discretización para pintar la analítica
 +
 +
% Condiciones iniciales
 +
u = 10 * ones(Nx, 1);
 +
 +
% Definimos la solucón estacionaria
 +
 +
estacionaria = @(x)9*x +1;
 +
 +
% Condiciones de frontera
 +
u(1) = 1;  % Extremo izquierdo
 +
u(end) = 10; % Extremo derecho
 +
 +
% Coeficiente del método
 +
r = alpha * dt / dx^2;
 +
if r > 0.5
 +
    error('El método no es estable, reduce dt o aumenta dx.');
 +
end
 +
 +
% Preparar la grabación de un GIF
 +
filename = 'heat_equation_comparacion_2.gif';
 +
frame_delay = 0.1;
 +
 +
% Evolución temporal
 +
for n = 1:Nt
 +
    u_new = u; % Copia de la solución actual
 +
    for i = 2:Nx-1
 +
        u_new(i) = u(i) + r * (u(i+1) - 2*u(i) + u(i-1));
 +
    end
 +
    u = u_new;
 +
   
 +
    % Solución analítica
 +
    t = n * dt;
 +
    w = zeros(size(xx));
 +
    for k = 1:100 % Consideramos los primeros 100 términos
 +
        w = w + (18 / (k * pi)) * sin(k * pi * xx) * exp(-k^2 * pi^2 * t);
 +
%        w = w + (-(18*(sin(pi*k) - pi*k))/(k^2*pi^2))*sin(k * pi * x) * exp(-k^2 * pi^2 * t);
 +
 +
    end
 +
   
 +
    % Guardar frames en el GIF
 +
    if mod(n, 10) == 0
 +
        plot(x, u, 'b-', 'LineWidth', 2);
 +
        hold on;
 +
        plot(xx, w + 9*xx + 1, 'r--', 'LineWidth', 2); % +1 para respetar la condición de frontera izquierda
 +
        plot(xx, estacionaria(xx), 'g-', 'LineWidth', 2); % Solución estacionaria
 +
        hold off;
 +
        xlabel('Posición'); ylabel('Temperatura');
 +
        title(sprintf('Tiempo: %.3f s', n*dt));
 +
        legend('Diferencias Finitas', 'Solución Analítica', 'Solución estacionaria');
 +
        drawnow;
 +
        frame = getframe(gcf);
 +
        img = frame2im(frame);
 +
        [A, map] = rgb2ind(img, 256);
 +
        if n == 10
 +
            imwrite(A, map, filename, 'gif', 'LoopCount', inf, 'DelayTime', frame_delay);
 +
        else
 +
            imwrite(A, map, filename, 'gif', 'WriteMode', 'append', 'DelayTime', frame_delay);
 +
        end
 +
    end
 +
end
 +
 +
% Gráfica final
 +
figure;
 +
plot(x, u, 'b-', 'LineWidth', 2);
 +
hold on;
 +
plot(xx, w + 9*xx+ 1, 'r--', 'LineWidth', 2);
 +
plot(xx, estacionaria(xx), 'g-', 'LineWidth', 2);
 +
hold off;
 +
xlabel('Posición'); ylabel('Temperatura');
 +
title('Comparación final de temperatura en la varilla');
 +
legend('Diferencias Finitas', 'Solución Analítica', 'Solución estacionaria');
 +
 +
}}
 +
 +
Observamos que la solución numérica aproxima mejor según avanza el tiempo. Analizaremos este fenómeno mediante los errores  en norma \( ||\cdot||_\infty \).
 +
 +
[[Archivo:Primeraimagen_ACIRV.jpg|600px|thumb|right|Error absoluto]]
 +
 +
{{matlab|codigo=
 +
% Resolución de la ecuación del calor 1D por diferencias finitas y comparación con solución analítica
 +
clear all; close all
 +
% Parámetros
 +
L = 1;        % Longitud de la varilla
 +
Nx = 4;      % Número de puntos en el espacio
 +
dx = L / (Nx-1);
 +
alpha = 1; % Difusividad térmica
 +
dt = 0.0005;  % Paso de tiempo
 +
Nt = 1000;    % Número de pasos de tiempo
 +
err = zeros(length(Nt));
 +
 +
% Discretización espacial
 +
x = linspace(0, L, Nx); % Discretización para la numérica
 +
xx = linspace(0, L, Nx); % Usamos la misma malla para comparación
 +
tt = linspace(0,Nt*dt,Nt);
 +
 +
% Condiciones iniciales
 +
u = 10 * ones(Nx, 1);
 +
 +
% Definimos la solución estacionaria
 +
estacionaria = @(x) 9*x + 1;
 +
 +
% Condiciones de frontera
 +
u(1) = 1;  % Extremo izquierdo
 +
u(end) = 10; % Extremo derecho
 +
 +
% Coeficiente del método
 +
r = alpha * dt / dx^2;
 +
if r > 0.5
 +
    error('El método no es estable, reduce dt o aumenta dx.');
 +
end
 +
 +
% Evolución temporal
 +
for n = 1:Nt
 +
    u_new = u;
 +
    for i = 2:Nx-1
 +
        u_new(i) = u(i) + r * (u(i+1) - 2*u(i) + u(i-1));
 +
    end
 +
    u = u_new;
 +
    w = zeros(size(xx));
 +
    t = n * dt;
 +
    for k = 1:100 % 100 términos en la serie
 +
        w = w + (18 / (k * pi)) * sin(k * pi * xx) * exp(-k^2 * pi^2 * t);
 +
    end
 +
    sol_analitica = w + 9*xx + 1;
 +
    err(n) = max(abs(sol_analitica-u'));
 +
end
 +
 +
% Cálculo de la solución analítica en el último tiempo
 +
w = zeros(size(xx));
 +
t = Nt * dt;
 +
for k = 1:100 % 100 términos en la serie
 +
    w = w + (18 / (k * pi)) * sin(k * pi * xx) * exp(-k^2 * pi^2 * t);
 +
end
 +
sol_analitica = w + 9*xx + 1;
 +
 +
% Cálculo del error
 +
error_abs = abs(u - sol_analitica');
 +
error_rel = error_abs ./ abs(sol_analitica');
 +
 +
% Gráficas
 +
figure;
 +
% subplot(2,1,1);
 +
% plot(x, u, 'b-', 'LineWidth', 2);
 +
% hold on;
 +
% plot(xx, sol_analitica, 'r--', 'LineWidth', 2);
 +
% plot(xx, estacionaria(xx), 'g-', 'LineWidth', 2);
 +
% hold off;
 +
% xlabel('Posición'); ylabel('Temperatura');
 +
% title('Comparación final de temperatura en la varilla');
 +
% legend('Diferencias Finitas', 'Solución Analítica', 'Solución estacionaria');
 +
%
 +
% % Gráfica del error
 +
% subplot(2,1,2);
 +
plot(tt, err, 'bl-', 'LineWidth', 2);
 +
xlabel('Posición'); ylabel('Error absoluto');
 +
plot(tt, log10(err), 'bl-', 'LineWidth', 2);
 +
xlabel('Posición'); ylabel('Error absoluto logaritmico');
 +
 +
title(sprintf('Error absoluto entre solución numérica y analítica para N=%d', Nx-1));
 +
grid on;
 +
 +
}}
 +
 +
Observamos que este error, sin ser muy significativo, es del orden de \(10^{-2}\) hablando en términos relativos. Destaca que comience siendo muy pequeño, esto se abordará al final de la sección. Veamos como mejorarlo al aumentar la malla espacial hasta \(N = 19 \) (20 puntos).
 +
 +
[[Archivo:Heat equation comparacion1.gif|600px|thumb|right|Comparación solución analítica y numérica para 20 puntos]]
 +
 +
{{matlab|codigo=
 +
% Resolución de la ecuación del calor 1D por diferencias finitas y comparación con solución analítica
 +
 +
% Parámetros
 +
L = 1;        % Longitud de la varilla
 +
Nx = 20;      % Número de puntos en el espacio
 +
dx = L / (Nx-1);
 +
alpha = 1; % Difusividad térmica
 +
dt = 0.0005;  % Paso de tiempo
 +
Nt = 1000;    % Número de pasos de tiempo
 +
 +
% Discretización espacial
 +
x = linspace(0, L, Nx); %Discretización para hacer los cálculos con la numérica
 +
xx = linspace(0,L,50); %Discretización para pintar la analítica
 +
 +
% Condiciones iniciales
 +
u = 10 * ones(Nx, 1);
 +
 +
% Definimos la solucón estacionaria
 +
 +
estacionaria = @(x)9*x +1;
 +
 +
% Condiciones de frontera
 +
u(1) = 1;  % Extremo izquierdo
 +
u(end) = 10; % Extremo derecho
 +
 +
% Coeficiente del método
 +
r = alpha * dt / dx^2;
 +
if r > 0.5
 +
    error('El método no es estable, reduce dt o aumenta dx.');
 +
end
 +
 +
% Preparar la grabación de un GIF
 +
filename = 'heat_equation_comparacion1.gif';
 +
frame_delay = 0.1;
 +
 +
% Evolución temporal
 +
for n = 1:Nt
 +
    u_new = u; % Copia de la solución actual
 +
    for i = 2:Nx-1
 +
        u_new(i) = u(i) + r * (u(i+1) - 2*u(i) + u(i-1));
 +
    end
 +
    u = u_new;
 +
   
 +
    % Solución analítica
 +
    t = n * dt;
 +
    w = zeros(size(xx));
 +
    for k = 1:100 % Consideramos los primeros 100 términos
 +
        w = w + (18 / (k * pi)) * sin(k * pi * xx) * exp(-k^2 * pi^2 * t);
 +
%        w = w + (-(18*(sin(pi*k) - pi*k))/(k^2*pi^2))*sin(k * pi * x) * exp(-k^2 * pi^2 * t);
 +
 +
    end
 +
   
 +
    % Guardar frames en el GIF
 +
    if mod(n, 10) == 0
 +
        plot(x, u, 'b-', 'LineWidth', 2);
 +
        hold on;
 +
        plot(xx, w + 9*xx + 1, 'r--', 'LineWidth', 2); % +1 para respetar la condición de frontera izquierda
 +
        plot(xx, estacionaria(xx), 'g-', 'LineWidth', 2); % Solución estacionaria
 +
        hold off;
 +
        xlabel('Posición'); ylabel('Temperatura');
 +
        title(sprintf('Tiempo: %.3f s', n*dt));
 +
        legend('Diferencias Finitas', 'Solución Analítica', 'Solución estacionaria');
 +
        drawnow;
 +
        frame = getframe(gcf);
 +
        img = frame2im(frame);
 +
        [A, map] = rgb2ind(img, 256);
 +
        if n == 10
 +
            imwrite(A, map, filename, 'gif', 'LoopCount', inf, 'DelayTime', frame_delay);
 +
        else
 +
            imwrite(A, map, filename, 'gif', 'WriteMode', 'append', 'DelayTime', frame_delay);
 +
        end
 +
    end
 +
end
 +
 +
% Gráfica final
 +
figure;
 +
plot(x, u, 'b-', 'LineWidth', 2);
 +
hold on;
 +
plot(xx, w + 9*xx+ 1, 'r--', 'LineWidth', 2);
 +
plot(xx, estacionaria(xx), 'g-', 'LineWidth', 2);
 +
hold off;
 +
xlabel('Posición'); ylabel('Temperatura');
 +
title('Comparación final de temperatura en la varilla');
 +
legend('Diferencias Finitas', 'Solución Analítica', 'Solución estacionaria');
 +
 +
}}
 +
 +
Se puede ver que se ajusta muy bien a la solución analítica. Véase en los errores, que precisan de una transformación logarítmica debido a la escala; graficaremos también la escala original para comparar con el anterior.
 +
 +
[[Archivo:Segundaimagen_ACIRV.jpg|600px|thumb|right|Error absoluto]]
 +
[[Archivo:logaritmica_3_ACIRV.jpg|600px|thumb|right|Error absoluto logarítmico]]
 +
 +
{{matlab|codigo=
 +
% Resolución de la ecuación del calor 1D por diferencias finitas y comparación con solución analítica
 +
clear all; close all
 +
% Parámetros
 +
L = 1;        % Longitud de la varilla
 +
Nx = 20;      % Número de puntos en el espacio
 +
dx = L / (Nx-1);
 +
alpha = 1; % Difusividad térmica
 +
dt = 0.0005;  % Paso de tiempo
 +
Nt = 1000;    % Número de pasos de tiempo
 +
err = zeros(length(Nt));
 +
 +
% Discretización espacial
 +
x = linspace(0, L, Nx); % Discretización para la numérica
 +
xx = linspace(0, L, Nx); % Usamos la misma malla para comparación
 +
tt = linspace(0,Nt*dt,Nt);
 +
 +
% Condiciones iniciales
 +
u = 10 * ones(Nx, 1);
 +
 +
% Definimos la solución estacionaria
 +
estacionaria = @(x) 9*x + 1;
 +
 +
% Condiciones de frontera
 +
u(1) = 1;  % Extremo izquierdo
 +
u(end) = 10; % Extremo derecho
 +
 +
% Coeficiente del método
 +
r = alpha * dt / dx^2;
 +
if r > 0.5
 +
    error('El método no es estable, reduce dt o aumenta dx.');
 +
end
 +
 +
% Evolución temporal
 +
for n = 1:Nt
 +
    u_new = u;
 +
    for i = 2:Nx-1
 +
        u_new(i) = u(i) + r * (u(i+1) - 2*u(i) + u(i-1));
 +
    end
 +
    u = u_new;
 +
    w = zeros(size(xx));
 +
    t = n * dt;
 +
    for k = 1:100 % 100 términos en la serie
 +
        w = w + (18 / (k * pi)) * sin(k * pi * xx) * exp(-k^2 * pi^2 * t);
 +
    end
 +
    sol_analitica = w + 9*xx + 1;
 +
    err(n) = max(abs(sol_analitica-u'));
 +
end
 +
 +
% Cálculo de la solución analítica en el último tiempo
 +
w = zeros(size(xx));
 +
t = Nt * dt;
 +
for k = 1:100 % 100 términos en la serie
 +
    w = w + (18 / (k * pi)) * sin(k * pi * xx) * exp(-k^2 * pi^2 * t);
 +
end
 +
sol_analitica = w + 9*xx + 1;
 +
 +
% Cálculo del error
 +
error_abs = abs(u - sol_analitica');
 +
error_rel = error_abs ./ abs(sol_analitica');
 +
 +
% Gráficas
 +
figure;
 +
% subplot(2,1,1);
 +
% plot(x, u, 'b-', 'LineWidth', 2);
 +
% hold on;
 +
% plot(xx, sol_analitica, 'r--', 'LineWidth', 2);
 +
% plot(xx, estacionaria(xx), 'g-', 'LineWidth', 2);
 +
% hold off;
 +
% xlabel('Posición'); ylabel('Temperatura');
 +
% title('Comparación final de temperatura en la varilla');
 +
% legend('Diferencias Finitas', 'Solución Analítica', 'Solución estacionaria');
 +
%
 +
% % Gráfica del error
 +
% subplot(2,1,2);
 +
plot(tt, err, 'bl-', 'LineWidth', 2);
 +
xlabel('Posición'); ylabel('Error absoluto');
 +
plot(tt, log10(err), 'bl-', 'LineWidth', 2);
 +
xlabel('Posición'); ylabel('Error absoluto logaritmico');
  
 +
title(sprintf('Error absoluto entre solución numérica y analítica para N=%d', Nx-1));
 +
grid on;
  
 +
}}
  
 +
Nótese que el error es muy pequeño, con orden de \(10^{-4}\). Destaca que para \(t\) cercanas a \(0\) el error sea alto, contrariamente a lo ocurrido anteriormente. La explicación viene de la aproximación de la serie de Fourier al término 100. Veámoslo gráficamente:
  
 +
[[Archivo:ACIRV tiempo inicial 1.jpg|500px|thumb|right|4 puntos en t = 0 con solución analítica]]
 +
[[Archivo:ACIRV tiempo inicial 2.jpg|500px|thumb|right|20 puntos en t = 0 con solución analítica]]
 +
{{matlab|codigo=
  
 +
L = 1;        % Longitud de la varilla
 +
Nx = 20;      % Número de puntos en el espacio
 +
dx = L / (Nx-1);
 +
alpha = 1; % Difusividad térmica
 +
dt = 0.0005;  % Paso de tiempo
 +
Nt = 1000;    % Número de pasos de tiempo
  
 +
% Discretización espacial
 +
x = linspace(0, L, Nx); %Discretización para hacer los cálculos con la numérica
 +
xx = linspace(0,L,50); %Discretización para pintar la analítica
  
 +
% Condiciones iniciales
 +
u = 10 * ones(Nx, 1);
  
 +
% Definimos la solucón estacionaria
  
 +
estacionaria = @(x)9*x +1;
  
 +
% Condiciones de frontera
 +
u(1) = 1;  % Extremo izquierdo
 +
u(end) = 10; % Extremo derecho
  
 +
% Coeficiente del método
 +
r = alpha * dt / dx^2;
 +
if r > 0.5
 +
    error('El método no es estable, reduce dt o aumenta dx.');
 +
end
  
  
 +
% Evolución temporal
 +
% Solución analítica
 +
t = 0;
 +
w = zeros(size(xx));
 +
for k = 1:100 % Consideramos los primeros 100 términos
 +
    w = w + (18 / (k * pi)) * sin(k * pi * xx) * exp(-k^2 * pi^2 * t);
 +
%        w = w + (-(18*(sin(pi*k) - pi*k))/(k^2*pi^2))*sin(k * pi * x) * exp(-k^2 * pi^2 * t);
  
 +
end
  
 +
% Guardar frames en el GIF
 +
plot(x, u, 'b-', 'LineWidth', 2);
 +
hold on;
 +
plot(x, u, 'bo', 'LineWidth', 3);
  
 +
plot(xx, w + 9*xx + 1, 'r--', 'LineWidth', 2); % +1 para respetar la condición de frontera izquierda
 +
plot(xx, estacionaria(xx), 'g-', 'LineWidth', 2); % Solución estacionaria
 +
hold off;
 +
xlabel('Posición'); ylabel('Temperatura');
 +
title(sprintf('Tiempo: %.3f s', 0));
 +
legend('Diferencias Finitas','Puntos sol numérica', 'Solución Analítica', 'Solución estacionaria');
 +
drawnow;
  
 +
}}
  
 +
La función analítica debería tener temperatura \(= 10\) para \(x \in (0,1]\), pero la aproximación para \(t = 0\) es mala en ambos casos. En el primero, todos los puntos quedan cercanos a la función, pero la aproximación dista mucho de la función analítica. En el segundo, el segundo punto se aleja de la función analítica por la aproximación pero los siguientes se acercan rápidamente; para infinitos términos de la serie esto se corrige.
  
 
[[Categoría:EDP]]
 
[[Categoría:EDP]]
 
[[Categoría:EDP24/25]]
 
[[Categoría:EDP24/25]]

Revisión actual del 00:18 20 mar 2025

Trabajo realizado por estudiantes
Título Ecuación del calor (Grupo ACIRV).
Asignatura EDP
Curso 2024-25
Autores Ángela Sotelo Fernández, Carmen Doñoro Molina, Inés Torres Gómez, Rubén Gutiérrez Hernández, Violeta Luján Barrios.
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura


1 Introducción

En esta página ilustraremos el uso del método de diferencias finitas para resolver la ecuación del calor en una dimensión.

Dicho método sirve para resolver de manera aproximada ecuaciones diferenciales parciales y ordinarias de manera numérica. Funciona sustituyendo las derivadas por cocientes de diferencias que vienen dados por valores de la función en puntos discretos de una malla.


2 Modelización de la ecuación del calor con una dimensión

La ecuación a resolver es

[math] u_t = \alpha u_{xx}, \quad x \in [0,1], \ t \gt 0,\ \alpha = 1. [/math]

donde \(u(x,t)\) representa la temperatura en el punto \(x\) y en el instante \(t\), y \(\alpha\) es una constante que viene dada por la conductividad térmica del material.

Se considera una varilla metálica en el intervalo \([0,1]\) y aislada en su superficie lateral tal que la conducción de calor solo se produce en la dirección longitudinal. La temperatura inicial de la varilla es \(10^{\circ}C\). En los extremos derecho e izquierdo la temperatura se mantiene a \(10^{\circ}C\) y \(1^{\circ}C\) respectivamente.

El sistema que modeliza el comportamiento de la temperatura, representada por la función \(u(x,t)\), es el siguiente:

[math] \left \{ \begin{array}{llll} u_t - u_{xx} & = 0 & x \in [0,1], & t\gt0 \\ u(0,t)& = 1 & t\gt0\\ u(1,t) & = 10 & t\gt0\\ u(x,0) &= 10 & x \in [0,1] \end{array} \right . \tag{1} [/math]

3 Resolución analítica del problema

Primeramente, resolveremos el problema analíticamente para después comparar su solución con la obtenida numéricamente.

Puesto que el sistema obtenido \((1)\) no es homogéneo, debemos encontrar primero la solución estacionaria. Es decir, suponemos que \(t \rightarrow \infty\) (la solución ya no varía en el tiempo).

La solución estacionaria obtenida es \(v(x) = 9x +1\) . Al representarla gráficamente en 3D obtenemos:

Solución estacionaria
%Creamos una matriz que representa una malla de puntos (tiempo,
% espacio) y en cada una de las columnas introducimos el valor de 9x+1 para todos los tiempos.

% Matriz
evaluaciones = zeros(100,100);

% Límite temporal y vectores 
T = 1;
t = linspace(0,T,100);
x = linspace(0,1,100);

% Rellenamos la matriz 
for j = 1:100
    evaluaciones(:,j) = (9*x(j) + 1) * ones(100,1);
end

% Representación gráfica
surf(t,x,evaluaciones')
title('Solución estacionaria: v(x) = 9x + 1')
xlabel('Tiempo')
ylabel('Espacio')
zlabel('Temperatura (°C)')


Una vez calculada la solución estacionaria, homogeneizamos el sistema definiendo la ecuación \(w(x,t) = u(x,t) - v(x)\).

[math] \left \{ \begin{array}{llll} w_t - w_{xx} & = 0 & x \in [0,1], & t\gt0 \\ w(0,t)& = 0 & t\gt0\\ w(1,t) & = 0 & t\gt0\\ w(x,0) &= 9(1-x) & x \in [0,1] \end{array} \right . [/math]

Resolviendo mediante el método de separación de variables y por superposición de soluciones, obtenemos la siguiente solución:

[math] w(x,t) = \sum^{\infty}_{k=1} \frac{18}{k\pi}sen(k\pi x)e^{-k^2 \pi^2 t} . [/math]

Véase la evolución de la temperatura en función del tiempo en cada punto de la barra. Representamos la temperatura respecto de la posición, variando el tiempo hasta \(t = 0.5 s\). Podemos representarla para diferentes términos, veamos que utilizando 100 se aproxima suficiente a tomar infinitos términos.


Solución analítica
Solución analítica 2 términos
Solución analítica 10 términos
% Resolución de la ecuación del calor 1D por diferencias finitas y comparación con solución analítica
clear all; close all;
% Parámetros
L = 1;        % Longitud de la varilla
Nx = 20;      % Número de puntos en el espacio
dx = L / (Nx-1);
alpha = 1; % Difusividad térmica
dt = 0.0005;  % Paso de tiempo
Nt = 1000;     % Número de pasos de tiempo

% Discretización espacial
x = linspace(0, L, Nx);

% Condiciones iniciales
u = 10 * ones(Nx, 1); 

% Condiciones de frontera
u(1) = 1;   % Extremo izquierdo
u(end) = 10; % Extremo derecho

% Definimos la solucón estacionaria

estacionaria = @(x)9*x +1;

% Coeficiente del método
r = alpha * dt / dx^2;
if r > 0.5
    error('El método no es estable, reduce dt o aumenta dx.');
end

% Preparar la grabación de un GIF
filename = 'heat_equation_analitico.gif';
frame_delay = 0.1;

% Evolución temporal
for n = 1:Nt
    % Solución analítica
    t = n * dt;
    w = zeros(size(x));
    for k = 1:100 % Consideramos los primeros 100 términos
%         w = w + (18 / (k * pi)) * sin(k * pi * x) * exp(-k^2 * pi^2 * t);
        w = w + (-(18*(sin(pi*k) - pi*k))/(k^2*pi^2))*sin(k * pi * x) * exp(-k^2 * pi^2 * t);

    end
    
    % Guardar frames en el GIF
    if mod(n, 10) == 0
        plot(x, estacionaria(x), 'g-', 'LineWidth', 2);
        hold on;
        plot(x, w + 9*x+ 1, 'r--', 'LineWidth', 2); % +1 para respetar la condición de frontera izquierda
        hold off;
        xlabel('Posición'); ylabel('Temperatura');
        title(sprintf('Tiempo: %.3f s', n*dt));
        legend('Solución estacionaria','Solución Analítica');
        drawnow;
        frame = getframe(gcf);
        img = frame2im(frame);
        [A, map] = rgb2ind(img, 256);
        if n == 10
            imwrite(A, map, filename, 'gif', 'LoopCount', inf, 'DelayTime', frame_delay);
        else
            imwrite(A, map, filename, 'gif', 'WriteMode', 'append', 'DelayTime', frame_delay);
        end
    end
end

% Gráfica final
figure;
plot(x, estacionaria(x), 'g-', 'LineWidth', 2);
hold on;
plot(x, w + 9*x + 1, 'r--', 'LineWidth', 2);
hold off;
xlabel('Posición'); ylabel('Temperatura');
title('Comparación final de temperatura en la varilla');
legend('Solución estacionaria', 'Solución Analítica');


Vemos que según aumenta el tiempo, se ajusta a la solución estacionaria.

4 Resolución numérica de la ecuación del calor

Falta obtener la resolución numérica del sistema planteado anteriormente \((1)\).

Así, discretizaremos el espacio en una malla de puntos espaciales \( x_i \) y tiempos \( t_n \). Estos deben cumplir \(x_0 = 0, \ x_N = 1\) de acuerdo con los intervalos de definición de la EDP. Definimos \(\Delta t = t_{n+1} - t_n\) y \(\Delta x = x_{i+1} - x_i\).

Pasamos a discretizar la ecuación según el método de diferencias finitas explícito. De \(u_t = \alpha u_{xx}\), pasamos a

[math] \frac{u_{i}^{n+1} - u_{i}^{n}}{\Delta t} = \alpha \frac{u_{i+1}^{n} - 2u_{i}^{n} + u_{i-1}^{n}}{\Delta x^2}. [/math]

Despejando \( u_{i}^{n+1} \), obtenemos la ecuación de actualización

[math] u_{i}^{n+1} = u_{i}^{n} + r \left( u_{i+1}^{n} - 2 u_{i}^{n} + u_{i-1}^{n}\right) . [/math]

Definimos el número de Fourier \(r = \frac{\alpha \Delta t}{\Delta x^2} \), que mide la relación entre la difusión térmica, el paso de tiempo y el desplazamiento. En nuestro caso tenemos \(\alpha = 1\).


Por último, inicializamos el problema para las condiciones de contorno. Definimos \( u_0^n = 1 , u_N^n = 10 \). Y las condiciones iniciales \(u_i^0 = 10, \ \forall i = 1,...,N.\)

Vamos a representar la solución numérica para diferentes mallados


Solución numérica para N = 3 (4 puntos y \(\Delta x = 0.33\)), \(\Delta t = 0.0005\)
Solución numérica para N = 20 (21 puntos y \(\Delta x = 0.05\)), \(\Delta t = 0.0005\)
% Resolución de la ecuación del calor 1D por diferencias finitas

% Parámetros
L = 1;        % Longitud de la varilla
Nx = 21;      % Número de puntos en el espacio
dx = L / (Nx-1);
D = 1; % Difusividad térmica
dt = 0.0005;  % Paso de tiempo
Nt = 1000;     % Número de pasos de tiempo

% Discretización espacial
x = linspace(0, L, Nx);

% Condiciones iniciales
u = 10 * ones(Nx, 1); 

% Condiciones de frontera
u(1) = 1;   % Extremo izquierdo
u(end) = 10; % Extremo derecho

%Solución estacionaria
estacionaria = @(x)9*x +1;


% Coeficiente del método
r = D * dt / dx^2;
if r > 0.5
    error('El método no es estable, reduce dt o aumenta dx.');
end

% Preparar la grabación de un GIF
filename = 'heat_equation_num_20.gif';
frame_delay = 0.1;

% Evolución temporal
for n = 1:Nt
    u_new = u; % Copia de la solución actual
    for i = 2:Nx-1
        u_new(i) = u(i) + r * (u(i+1) - 2*u(i) + u(i-1));
    end
    u = u_new;
    
    % Guardar frames en el GIF
    if mod(n, 10) == 0
        plot(x, u, 'b-', 'LineWidth', 2);
        hold on;
        plot(xx, estacionaria(xx), 'g-', 'LineWidth', 2); % Solución estacionaria
        hold off;        xlabel('Posición'); ylabel('Temperatura');
        title(sprintf('Tiempo: %.3f s', n*dt));
        legend('Diferencias Finitas', 'Solución estacionaria');
        drawnow;
        frame = getframe(gcf);
        img = frame2im(frame);
        [A, map] = rgb2ind(img, 256);
        if n == 10
            imwrite(A, map, filename, 'gif', 'LoopCount', inf, 'DelayTime', frame_delay);
        else
            imwrite(A, map, filename, 'gif', 'WriteMode', 'append', 'DelayTime', frame_delay);
        end
    end
end

% Gráfica final
figure;
plot(x, u, 'b-', 'LineWidth', 2);
xlabel('Posición'); ylabel('Temperatura');

title('Distribución final de temperatura en la varilla');


5 Análisis de estabilidad

Es importante asegurarse de que el método es estable. Aplicamos el análisis de von Neumann, en el que asumimos soluciones de la forma:

[math] u_{i}^{n} = G^n e^{ikx_i} , [/math]

donde \( G \) es el factor de amplificación. Reemplazando esta forma en la ecuación de actualización y tras un análisis algebraico, se obtiene:

[math] G = 1 - 4r \sin^2 \left( \frac{k \Delta x}{2} \right) . [/math]

Para que el método sea estable, es necesario que \( |G| \leq 1 \). Como el término \( \sin^2 (\cdot) \) está acotado entre 0 y 1, se deduce la condición de estabilidad:

[math] 1 - 4r \leq 1 \Rightarrow r \leq \frac{1}{2} . [/math]

Si \( r > \frac{1}{2} \), entonces \( |G| > 1 \), haciendo que la solución sea inestable.

6 Comparación de resultados y análisis de errores

Seguidamente, compararemos la solución analítica con la numérica para verificar la precisión de esta última en función de la malla empleada.

La siguiente representación recoge cómo varía la temperatura según qué solución en función del tiempo en cada punto.

Comparación solución analítica y numérica para 4 puntos
% Resolución de la ecuación del calor 1D por diferencias finitas y comparación con solución analítica

% Parámetros
L = 1;        % Longitud de la varilla
Nx = 4;      % Número de puntos en el espacio
dx = L / (Nx-1);
alpha = 1; % Difusividad térmica
dt = 0.0005;  % Paso de tiempo
Nt = 1000;     % Número de pasos de tiempo

% Discretización espacial
x = linspace(0, L, Nx); %Discretización para hacer los cálculos con la numérica
xx = linspace(0,L,50); %Discretización para pintar la analítica

% Condiciones iniciales
u = 10 * ones(Nx, 1); 

% Definimos la solucón estacionaria

estacionaria = @(x)9*x +1;

% Condiciones de frontera
u(1) = 1;   % Extremo izquierdo
u(end) = 10; % Extremo derecho

% Coeficiente del método
r = alpha * dt / dx^2;
if r > 0.5
    error('El método no es estable, reduce dt o aumenta dx.');
end

% Preparar la grabación de un GIF
filename = 'heat_equation_comparacion_2.gif';
frame_delay = 0.1;

% Evolución temporal
for n = 1:Nt
    u_new = u; % Copia de la solución actual
    for i = 2:Nx-1
        u_new(i) = u(i) + r * (u(i+1) - 2*u(i) + u(i-1));
    end
    u = u_new;
    
    % Solución analítica
    t = n * dt;
    w = zeros(size(xx));
    for k = 1:100 % Consideramos los primeros 100 términos
        w = w + (18 / (k * pi)) * sin(k * pi * xx) * exp(-k^2 * pi^2 * t);
%         w = w + (-(18*(sin(pi*k) - pi*k))/(k^2*pi^2))*sin(k * pi * x) * exp(-k^2 * pi^2 * t);

    end
    
    % Guardar frames en el GIF
    if mod(n, 10) == 0
        plot(x, u, 'b-', 'LineWidth', 2);
        hold on;
        plot(xx, w + 9*xx + 1, 'r--', 'LineWidth', 2); % +1 para respetar la condición de frontera izquierda
        plot(xx, estacionaria(xx), 'g-', 'LineWidth', 2); % Solución estacionaria
        hold off;
        xlabel('Posición'); ylabel('Temperatura');
        title(sprintf('Tiempo: %.3f s', n*dt));
        legend('Diferencias Finitas', 'Solución Analítica', 'Solución estacionaria');
        drawnow;
        frame = getframe(gcf);
        img = frame2im(frame);
        [A, map] = rgb2ind(img, 256);
        if n == 10
            imwrite(A, map, filename, 'gif', 'LoopCount', inf, 'DelayTime', frame_delay);
        else
            imwrite(A, map, filename, 'gif', 'WriteMode', 'append', 'DelayTime', frame_delay);
        end
    end
end

% Gráfica final
figure;
plot(x, u, 'b-', 'LineWidth', 2);
hold on;
plot(xx, w + 9*xx+ 1, 'r--', 'LineWidth', 2);
plot(xx, estacionaria(xx), 'g-', 'LineWidth', 2);
hold off;
xlabel('Posición'); ylabel('Temperatura');
title('Comparación final de temperatura en la varilla');
legend('Diferencias Finitas', 'Solución Analítica', 'Solución estacionaria');


Observamos que la solución numérica aproxima mejor según avanza el tiempo. Analizaremos este fenómeno mediante los errores en norma \( ||\cdot||_\infty \).

Error absoluto
% Resolución de la ecuación del calor 1D por diferencias finitas y comparación con solución analítica
clear all; close all
% Parámetros
L = 1;        % Longitud de la varilla
Nx = 4;      % Número de puntos en el espacio
dx = L / (Nx-1);
alpha = 1; % Difusividad térmica
dt = 0.0005;  % Paso de tiempo
Nt = 1000;     % Número de pasos de tiempo
err = zeros(length(Nt));

% Discretización espacial
x = linspace(0, L, Nx); % Discretización para la numérica
xx = linspace(0, L, Nx); % Usamos la misma malla para comparación
tt = linspace(0,Nt*dt,Nt);

% Condiciones iniciales
u = 10 * ones(Nx, 1); 

% Definimos la solución estacionaria
estacionaria = @(x) 9*x + 1;

% Condiciones de frontera
u(1) = 1;   % Extremo izquierdo
u(end) = 10; % Extremo derecho

% Coeficiente del método
r = alpha * dt / dx^2;
if r > 0.5
    error('El método no es estable, reduce dt o aumenta dx.');
end

% Evolución temporal
for n = 1:Nt
    u_new = u;
    for i = 2:Nx-1
        u_new(i) = u(i) + r * (u(i+1) - 2*u(i) + u(i-1));
    end
    u = u_new;
    w = zeros(size(xx));
    t = n * dt;
    for k = 1:100 % 100 términos en la serie
        w = w + (18 / (k * pi)) * sin(k * pi * xx) * exp(-k^2 * pi^2 * t);
    end
    sol_analitica = w + 9*xx + 1;
    err(n) = max(abs(sol_analitica-u'));
end

% Cálculo de la solución analítica en el último tiempo
w = zeros(size(xx));
t = Nt * dt;
for k = 1:100 % 100 términos en la serie
    w = w + (18 / (k * pi)) * sin(k * pi * xx) * exp(-k^2 * pi^2 * t);
end
sol_analitica = w + 9*xx + 1;

% Cálculo del error
error_abs = abs(u - sol_analitica');
error_rel = error_abs ./ abs(sol_analitica');

% Gráficas
figure;
% subplot(2,1,1);
% plot(x, u, 'b-', 'LineWidth', 2);
% hold on;
% plot(xx, sol_analitica, 'r--', 'LineWidth', 2);
% plot(xx, estacionaria(xx), 'g-', 'LineWidth', 2);
% hold off;
% xlabel('Posición'); ylabel('Temperatura');
% title('Comparación final de temperatura en la varilla');
% legend('Diferencias Finitas', 'Solución Analítica', 'Solución estacionaria');
% 
% % Gráfica del error
% subplot(2,1,2);
plot(tt, err, 'bl-', 'LineWidth', 2);
xlabel('Posición'); ylabel('Error absoluto');
plot(tt, log10(err), 'bl-', 'LineWidth', 2);
xlabel('Posición'); ylabel('Error absoluto logaritmico');

title(sprintf('Error absoluto entre solución numérica y analítica para N=%d', Nx-1));
grid on;


Observamos que este error, sin ser muy significativo, es del orden de \(10^{-2}\) hablando en términos relativos. Destaca que comience siendo muy pequeño, esto se abordará al final de la sección. Veamos como mejorarlo al aumentar la malla espacial hasta \(N = 19 \) (20 puntos).

Comparación solución analítica y numérica para 20 puntos
% Resolución de la ecuación del calor 1D por diferencias finitas y comparación con solución analítica

% Parámetros
L = 1;        % Longitud de la varilla
Nx = 20;      % Número de puntos en el espacio
dx = L / (Nx-1);
alpha = 1; % Difusividad térmica
dt = 0.0005;  % Paso de tiempo
Nt = 1000;     % Número de pasos de tiempo

% Discretización espacial
x = linspace(0, L, Nx); %Discretización para hacer los cálculos con la numérica
xx = linspace(0,L,50); %Discretización para pintar la analítica

% Condiciones iniciales
u = 10 * ones(Nx, 1); 

% Definimos la solucón estacionaria

estacionaria = @(x)9*x +1;

% Condiciones de frontera
u(1) = 1;   % Extremo izquierdo
u(end) = 10; % Extremo derecho

% Coeficiente del método
r = alpha * dt / dx^2;
if r > 0.5
    error('El método no es estable, reduce dt o aumenta dx.');
end

% Preparar la grabación de un GIF
filename = 'heat_equation_comparacion1.gif';
frame_delay = 0.1;

% Evolución temporal
for n = 1:Nt
    u_new = u; % Copia de la solución actual
    for i = 2:Nx-1
        u_new(i) = u(i) + r * (u(i+1) - 2*u(i) + u(i-1));
    end
    u = u_new;
    
    % Solución analítica
    t = n * dt;
    w = zeros(size(xx));
    for k = 1:100 % Consideramos los primeros 100 términos
        w = w + (18 / (k * pi)) * sin(k * pi * xx) * exp(-k^2 * pi^2 * t);
%         w = w + (-(18*(sin(pi*k) - pi*k))/(k^2*pi^2))*sin(k * pi * x) * exp(-k^2 * pi^2 * t);

    end
    
    % Guardar frames en el GIF
    if mod(n, 10) == 0
        plot(x, u, 'b-', 'LineWidth', 2);
        hold on;
        plot(xx, w + 9*xx + 1, 'r--', 'LineWidth', 2); % +1 para respetar la condición de frontera izquierda
        plot(xx, estacionaria(xx), 'g-', 'LineWidth', 2); % Solución estacionaria
        hold off;
        xlabel('Posición'); ylabel('Temperatura');
        title(sprintf('Tiempo: %.3f s', n*dt));
        legend('Diferencias Finitas', 'Solución Analítica', 'Solución estacionaria');
        drawnow;
        frame = getframe(gcf);
        img = frame2im(frame);
        [A, map] = rgb2ind(img, 256);
        if n == 10
            imwrite(A, map, filename, 'gif', 'LoopCount', inf, 'DelayTime', frame_delay);
        else
            imwrite(A, map, filename, 'gif', 'WriteMode', 'append', 'DelayTime', frame_delay);
        end
    end
end

% Gráfica final
figure;
plot(x, u, 'b-', 'LineWidth', 2);
hold on;
plot(xx, w + 9*xx+ 1, 'r--', 'LineWidth', 2);
plot(xx, estacionaria(xx), 'g-', 'LineWidth', 2);
hold off;
xlabel('Posición'); ylabel('Temperatura');
title('Comparación final de temperatura en la varilla');
legend('Diferencias Finitas', 'Solución Analítica', 'Solución estacionaria');


Se puede ver que se ajusta muy bien a la solución analítica. Véase en los errores, que precisan de una transformación logarítmica debido a la escala; graficaremos también la escala original para comparar con el anterior.

Error absoluto
Error absoluto logarítmico
% Resolución de la ecuación del calor 1D por diferencias finitas y comparación con solución analítica
clear all; close all
% Parámetros
L = 1;        % Longitud de la varilla
Nx = 20;      % Número de puntos en el espacio
dx = L / (Nx-1);
alpha = 1; % Difusividad térmica
dt = 0.0005;  % Paso de tiempo
Nt = 1000;     % Número de pasos de tiempo
err = zeros(length(Nt));

% Discretización espacial
x = linspace(0, L, Nx); % Discretización para la numérica
xx = linspace(0, L, Nx); % Usamos la misma malla para comparación
tt = linspace(0,Nt*dt,Nt);

% Condiciones iniciales
u = 10 * ones(Nx, 1); 

% Definimos la solución estacionaria
estacionaria = @(x) 9*x + 1;

% Condiciones de frontera
u(1) = 1;   % Extremo izquierdo
u(end) = 10; % Extremo derecho

% Coeficiente del método
r = alpha * dt / dx^2;
if r > 0.5
    error('El método no es estable, reduce dt o aumenta dx.');
end

% Evolución temporal
for n = 1:Nt
    u_new = u;
    for i = 2:Nx-1
        u_new(i) = u(i) + r * (u(i+1) - 2*u(i) + u(i-1));
    end
    u = u_new;
    w = zeros(size(xx));
    t = n * dt;
    for k = 1:100 % 100 términos en la serie
        w = w + (18 / (k * pi)) * sin(k * pi * xx) * exp(-k^2 * pi^2 * t);
    end
    sol_analitica = w + 9*xx + 1;
    err(n) = max(abs(sol_analitica-u'));
end

% Cálculo de la solución analítica en el último tiempo
w = zeros(size(xx));
t = Nt * dt;
for k = 1:100 % 100 términos en la serie
    w = w + (18 / (k * pi)) * sin(k * pi * xx) * exp(-k^2 * pi^2 * t);
end
sol_analitica = w + 9*xx + 1;

% Cálculo del error
error_abs = abs(u - sol_analitica');
error_rel = error_abs ./ abs(sol_analitica');

% Gráficas
figure;
% subplot(2,1,1);
% plot(x, u, 'b-', 'LineWidth', 2);
% hold on;
% plot(xx, sol_analitica, 'r--', 'LineWidth', 2);
% plot(xx, estacionaria(xx), 'g-', 'LineWidth', 2);
% hold off;
% xlabel('Posición'); ylabel('Temperatura');
% title('Comparación final de temperatura en la varilla');
% legend('Diferencias Finitas', 'Solución Analítica', 'Solución estacionaria');
% 
% % Gráfica del error
% subplot(2,1,2);
plot(tt, err, 'bl-', 'LineWidth', 2);
xlabel('Posición'); ylabel('Error absoluto');
plot(tt, log10(err), 'bl-', 'LineWidth', 2);
xlabel('Posición'); ylabel('Error absoluto logaritmico');

title(sprintf('Error absoluto entre solución numérica y analítica para N=%d', Nx-1));
grid on;


Nótese que el error es muy pequeño, con orden de \(10^{-4}\). Destaca que para \(t\) cercanas a \(0\) el error sea alto, contrariamente a lo ocurrido anteriormente. La explicación viene de la aproximación de la serie de Fourier al término 100. Veámoslo gráficamente:

4 puntos en t = 0 con solución analítica
20 puntos en t = 0 con solución analítica
L = 1;        % Longitud de la varilla
Nx = 20;      % Número de puntos en el espacio
dx = L / (Nx-1);
alpha = 1; % Difusividad térmica
dt = 0.0005;  % Paso de tiempo
Nt = 1000;     % Número de pasos de tiempo

% Discretización espacial
x = linspace(0, L, Nx); %Discretización para hacer los cálculos con la numérica
xx = linspace(0,L,50); %Discretización para pintar la analítica

% Condiciones iniciales
u = 10 * ones(Nx, 1); 

% Definimos la solucón estacionaria

estacionaria = @(x)9*x +1;

% Condiciones de frontera
u(1) = 1;   % Extremo izquierdo
u(end) = 10; % Extremo derecho

% Coeficiente del método
r = alpha * dt / dx^2;
if r > 0.5
    error('El método no es estable, reduce dt o aumenta dx.');
end


% Evolución temporal
% Solución analítica
t = 0;
w = zeros(size(xx));
for k = 1:100 % Consideramos los primeros 100 términos
    w = w + (18 / (k * pi)) * sin(k * pi * xx) * exp(-k^2 * pi^2 * t);
%         w = w + (-(18*(sin(pi*k) - pi*k))/(k^2*pi^2))*sin(k * pi * x) * exp(-k^2 * pi^2 * t);

end

% Guardar frames en el GIF
plot(x, u, 'b-', 'LineWidth', 2);
hold on;
plot(x, u, 'bo', 'LineWidth', 3);

plot(xx, w + 9*xx + 1, 'r--', 'LineWidth', 2); % +1 para respetar la condición de frontera izquierda
plot(xx, estacionaria(xx), 'g-', 'LineWidth', 2); % Solución estacionaria
hold off;
xlabel('Posición'); ylabel('Temperatura');
title(sprintf('Tiempo: %.3f s', 0));
legend('Diferencias Finitas','Puntos sol numérica', 'Solución Analítica', 'Solución estacionaria');
drawnow;


La función analítica debería tener temperatura \(= 10\) para \(x \in (0,1]\), pero la aproximación para \(t = 0\) es mala en ambos casos. En el primero, todos los puntos quedan cercanos a la función, pero la aproximación dista mucho de la función analítica. En el segundo, el segundo punto se aleja de la función analítica por la aproximación pero los siguientes se acercan rápidamente; para infinitos términos de la serie esto se corrige.

Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Ecuación_del_calor_(Grupo_ACIRV)&oldid=85128»