Ecuación del calor (Grupo ACIRV)

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:

Solución analítica

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.\)

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

Vamos a representar las soluciones obtenidas. Comenzamos con la analítica, calculada anteriormente. Representamos la temperatura respecto de la posición, variando el tiempo hasta \(t = 0.5 s\).

METER CÓDIGO

Solución analítica

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

A continuación, utilizamos también MATLAB para representar la solución numérica análogamente a la analítica. Podemos comparar ambas. Comenzamos con un mallado de 4 puntos en la posición ( \(N = 3 \) ). Para el temporal utilizamos \(1000 \) puntos y \(\Delta t = 0.0005 \), para que no introduzca error.

METER heat_equation_analitico_1 Y CÓDIGO

Observamos que aproxima mejor según se acerque a la estacionaria. Representamos los errores utilizando \( ||\cdot||_\infty \).

METER CÓDIGO

Error absoluto

Por último, mejoramos la malla espacial hasta \(N = 19 \) (20 puntos).

METER heat_equation_analitico_2 Y CÓDIGO

Se puede ver que ajusta en gran medida a la solución analítica. Observamos los errores.

METER CÓDIGO

Error absoluto

METER CÓDIGO

Error absoluto logarítmico