Diferencias finitas básicas

Discretización espacial y temporal de operadores diferenciales en mallas para simulación numérica.

1. Discretización de Dominios y Operadores Diferenciales

En física continua, las leyes fundamentales de conservación (masa, momento, energía, electromagnetismo) se expresan a través de Ecuaciones en Derivadas Parciales (EDPs). A diferencia de las EDOs, las EDPs describen campos escalares o vectoriales $u(\vec{x}, t)$ que varían simultáneamente en el espacio y en el tiempo.

El método de Diferencias Finitas (FDM) aproxima las derivadas continuas reemplazándolas por cocientes de diferencias evaluados sobre una malla discreta de puntos. Para un dominio unidimensional con longitud $L$ y tiempo $T$:

$x_i = i \Delta x, \quad i = 0, 1, \dots, N_x, \quad \Delta x = \frac{L}{N_x}$ $t_n = n \Delta t, \quad n = 0, 1, \dots, N_t, \quad \Delta t = \frac{T}{N_t}$

Se denota la aproximación numérica del campo en el nodo espacial $i$ y paso temporal $n$ como: $u_i^n \approx u(x_i, t_n)$

2. Fórmulas de Diferencias Finitas Derivadas de Taylor

A partir del desarrollo en serie de Taylor de $u(x \pm \Delta x)$:

$u(x + \Delta x) = u(x) + \Delta x \frac{\partial u}{\partial x} + \frac{\Delta x^2}{2} \frac{\partial^2 u}{\partial x^2} + \frac{\Delta x^3}{6} \frac{\partial^3 u}{\partial x^3} + \mathcal{O}(\Delta x^4)$ $u(x - \Delta x) = u(x) - \Delta x \frac{\partial u}{\partial x} + \frac{\Delta x^2}{2} \frac{\partial^2 u}{\partial x^2} - \frac{\Delta x^3}{6} \frac{\partial^3 u}{\partial x^3} + \mathcal{O}(\Delta x^4)$

Se deducen los operadores estándar:

1. Diferencia hacia adelante (Forward): $\left.\frac{\partial u}{\partial x}\right|_i = \frac{u_{i+1} - u_i}{\Delta x} + \mathcal{O}(\Delta x)$
2. Diferencia hacia atrás (Backward): $\left.\frac{\partial u}{\partial x}\right|_i = \frac{u_i - u_{i-1}}{\Delta x} + \mathcal{O}(\Delta x)$
3. Diferencia centrada (Central): $\left.\frac{\partial u}{\partial x}\right|_i = \frac{u_{i+1} - u_{i-1}}{2 \Delta x} + \mathcal{O}(\Delta x^2)$
4. Segunda derivada espacial centrada: $\left.\frac{\partial^2 u}{\partial x^2}\right|_i = \frac{u_{i+1} - 2u_i + u_{i-1}}{\Delta x^2} + \mathcal{O}(\Delta x^2)$

3. Ecuación de Difusión / Calor 1D y Esquema Explícito FTCS

La ecuación clásica de difusión térmica homogénea es: $\frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2}$ donde $\alpha$ es la difusividad térmica ($[\text{m}^2/\text{s}]$).

Aplicando diferencia hacia adelante en tiempo y centrada en espacio (FTCS: Forward-Time Central-Space): $\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}$

Despejando el valor futuro desconocido $u_i^{n+1}$: $u_i^{n+1} = u_i^n + r \left(u_{i+1}^n - 2u_i^n + u_{i-1}^n\right)$ donde el parámetro adimensional de difusión es: $r = \frac{\alpha \Delta t}{\Delta x^2}$

4. Análisis de Estabilidad de von Neumann

El análisis de von Neumann evalúa la evolución de un modo de Fourier arbitrario de la solución numérica: $u_i^n = G^n e^{i k x_i} = G^n e^{i k (i \Delta x)}$ donde $G$ es el factor de amplificación y $k$ es el número de onda espacial. Para que el esquema sea estable y los errores de redondeo no crezcan exponencialmente, se requiere: $|G| \le 1, \quad \forall k$

Sustituyendo el modo de Fourier en el esquema FTCS: $G^{n+1} e^{i k i \Delta x} = G^n e^{i k i \Delta x} \left[1 + r \left(e^{i k \Delta x} - 2 + e^{-i k \Delta x}\right)\right]$ $G = 1 + r \left(2 \cos(k \Delta x) - 2\right) = 1 - 4r \sin^2\left(\frac{k \Delta x}{2}\right)$

La condición $|G| \le 1$ impone: $-1 \le 1 - 4r \sin^2\left(\frac{k \Delta x}{2}\right) \le 1$ El límite más restrictivo ocurre cuando $\sin^2(k \Delta x / 2) = 1$: $-1 \le 1 - 4r \implies 4r \le 2 \implies r = \frac{\alpha \Delta t}{\Delta x^2} \le \frac{1}{2}$

Consecuencia Física Crucial:

Si se duplica la resolución espacial ($\Delta x \to \Delta x/2$), el paso temporal máximo permitido para mantener estabilidad disminuye por un factor de 4: $\Delta t \le \frac{\Delta x^2}{2\alpha}$

5. Ecuación de Onda 1D y Condición CFL

Para la propagación de ondas mecánicas o electromagnéticas en 1D: $\frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2}$ Discretizando ambas derivadas segundas mediante diferencias centradas: $\frac{u_i^{n+1} - 2u_i^n + u_i^{n-1}}{\Delta t^2} = c^2 \frac{u_{i+1}^n - 2u_i^n + u_{i-1}^n}{\Delta x^2}$

Despejando $u_i^{n+1}$: $u_i^{n+1} = 2u_i^n - u_i^{n-1} + C^2 \left(u_{i+1}^n - 2u_i^n + u_{i-1}^n\right)$ donde el Número de Courant adimensional es: $C = \frac{c \Delta t}{\Delta x}$

Condición de Courant-Friedrichs-Lewy (CFL):

El esquema es estable si y solo si: $C = \frac{c \Delta t}{\Delta x} \le 1$ Significado geométrico/físico: El dominio de dependencia numérico de la malla debe contener completamente el cono de luz o dominio de dependencia físico de la EDP hiperbólica. Si $C > 1$, la información física viaja más rápido que la malla, provocando inestabilidad numérica inmediata.

6. Esquemas Implícitos y Crank-Nicolson

Para eludir la severa restricción cuadrática $\Delta t \sim \Delta x^2$ en problemas de difusión, se emplean esquemas implícitos. El esquema de Crank-Nicolson promedia las derivadas espaciales en el tiempo actual $n$ y el tiempo futuro $n+1$:

$\frac{u_i^{n+1} - u_i^n}{\Delta t} = \frac{\alpha}{2} \left[\frac{u_{i+1}^{n+1} - 2u_i^{n+1} + u_{i-1}^{n+1}}{\Delta x^2} + \frac{u_{i+1}^n - 2u_i^n + u_{i-1}^n}{\Delta x^2}\right]$

• Precisión: Segundo orden en tiempo y en espacio $\mathcal{O}(\Delta t^2 + \Delta x^2)$.
• Estabilidad: Incondicionalmente estable para cualquier valor de $\Delta t$ ($|G| \le 1$ siempre).
• Costo: Requiere resolver un sistema tridiagonal lineal en cada paso temporal mediante el algoritmo de Thomas.

Procedimiento de Discretización y Resolución de EDPs por Diferencias Finitas

Metodología Numérica Sistemática

Para resolver computacionalmente un problema de valores iniciales y de contorno (IBVP) modelado por EDPs mediante diferencias finitas, se debe ejecutar el siguiente procedimiento:

Paso 1: Definición del Dominio y Malla Computacional

1. Especificar los límites espaciales $[0, L]$ y temporales $[0, T]$.
2. Escoger el número de subdivisiones espaciales $N_x$, determinando el paso $\Delta x = L / N_x$.
3. Localizar los nodos espaciales $x_i = i \Delta x$ para $i = 0, 1, \dots, N_x$.

Paso 2: Elección del Paso Temporal según Criterios de Estabilidad

• Para Ecuaciones Parabólicas (Calor/Difusión): $\Delta t \le \frac{\Delta x^2}{2 \alpha} \implies r = \frac{\alpha \Delta t}{\Delta x^2} \le 0.5$ Se recomienda fijar un margen de seguridad, por ejemplo $r = 0.4$.
• Para Ecuaciones Hiperbólicas (Ondas): $\Delta t \le \frac{\Delta x}{c} \implies C = \frac{c \Delta t}{\Delta x} \le 1.0$ Se recomienda fijar $C = 0.8 - 0.9$.

Paso 3: Inicialización de Condiciones Iniciales (t = 0)

1. Evaluar el perfil inicial continuo sobre la malla discreta: $u_i^0 = f(x_i), \quad \forall i = 0, 1, \dots, N_x$
2. En problemas de segundo orden temporal (ondas), donde se requiere $\dot{u}(x, 0) = g(x)$, aplicar un paso temporal ficticio hacia atrás para inicializar $u_i^1$: $u_i^1 = u_i^0 + \Delta t \, g(x_i) + \frac{C^2}{2} (u_{i+1}^0 - 2u_i^0 + u_{i-1}^0)$

Paso 4: Implementación de Condiciones de Frontera (BC)

En cada paso de tiempo $n$, aplicar rigurosamente las condiciones en los bordes $i = 0$ e $i = N_x$:

• Frontera de Dirichlet (Valor Fijo): $u_0^n = U_{\text{izq}}, \quad u_{N_x}^n = U_{\text{der}}$
• Frontera de Neumann (Flujo Fijo / Aislante $\partial u / \partial x = 0$): Utilizar un punto fantasma $u_{-1}^n$: $\frac{u_1^n - u_{-1}^n}{2 \Delta x} = 0 \implies u_{-1}^n = u_1^n$ Sustituyendo en el esquema interior: $u_0^{n+1} = u_0^n + 2r (u_1^n - u_0^n)$

Paso 5: Bucle Temporal y Actualización Espacial

1. Para cada paso $n = 0, 1, \dots, N_t - 1$:
   • Aplicar el stencil explícito en todos los nodos interiores $i = 1, 2, \dots, N_x - 1$.
   • Imponer las condiciones de frontera para el nuevo nivel $u^{n+1}$.
   • Intercambiar los buffers de memoria ($u^n \leftarrow u^{n+1}$).
2. Registrar diagnósticos físicos periódicos (energía elástica, flujo total de calor, masa conservada).

Ejemplos Resueltos: Difusión de Calor y Ondas 1D

Ejemplo 1: Conducción Térmica Transitoria en una Barra Unidimensional

Considérese una barra metálica de longitud $L = 1.0\text{ m}$ con difusividad térmica $\alpha = 1.0\text{ m}^2/\text{s}$. Las condiciones de frontera son extremos mantenidos a temperatura cero: $u(0, t) = 0, \quad u(1, t) = 0$ La distribución inicial de temperatura es: $u(x, 0) = 100 \sin(\pi x)$ La solución analítica exacta es: $u_{\text{exacta}}(x, t) = 100 e^{-\pi^2 \alpha t} \sin(\pi x)$

Discretización Numérica:

• Dividimos la barra en $N_x = 4$ intervalos iguales $\implies \Delta x = 0.25\text{ m}$.
• Los nodos espaciales son:
  • $x_0 = 0.0$, $x_1 = 0.25$, $x_2 = 0.50$, $x_3 = 0.75$, $x_4 = 1.0$.
• Temperaturas iniciales en $t = 0$:
  • $u_0^0 = 0$
  • $u_1^0 = 100 \sin(\pi/4) = 70.7107$
  • $u_2^0 = 100 \sin(\pi/2) = 100.0000$
  • $u_3^0 = 100 \sin(3\pi/4) = 70.7107$
  • $u_4^0 = 0$

Elegimos un parámetro de difusión estable $r = \frac{\alpha \Delta t}{\Delta x^2} = 0.40 \le 0.5$. Paso de tiempo: $\Delta t = r \frac{\Delta x^2}{\alpha} = 0.40 \times (0.25)^2 / 1.0 = 0.025\text{ s}$.

Cálculo del Primer Paso Temporal ($n = 0 \to n = 1$, $t = 0.025\text{ s}$):

Aplicamos el esquema FTCS: $u_i^1 = u_i^0 + 0.4 (u_{i+1}^0 - 2u_i^0 + u_{i-1}^0)$:

• Nodo 1 ($x_1 = 0.25$): $u_1^1 = 70.7107 + 0.4 [100.0000 - 2(70.7107) + 0.0] = 70.7107 + 0.4 [100.0 - 141.4214] = 70.7107 + 0.4(-41.4214) = 54.1421$

• Nodo 2 ($x_2 = 0.50$): $u_2^1 = 100.0000 + 0.4 [70.7107 - 2(100.0000) + 70.7107] = 100.0000 + 0.4 [141.4214 - 200.0] = 100.0000 + 0.4(-58.5786) = 76.5685$

• Nodo 3 ($x_3 = 0.75$): Por simetría física: $u_3^1 = u_1^1 = 54.1421$.

Comparación con la Solución Analítica en $t = 0.025\text{ s}$:

$u_{\text{exacta}}(0.5, 0.025) = 100 e^{-\pi^2 (1.0)(0.025)} \sin(\pi / 2) = 100 e^{-0.246740} = 78.1343$ El cálculo numérico con una malla extremadamente gruesa de solo 4 celdas arroja $76.5685$, con un error relativo de apenas el $2.0\%$.

Ejemplo 2: Propagación de un Pulso en una Cuerda Tensada (Ecuación de Onda)

Considérese una cuerda de longitud $L = 1\text{ m}$ con velocidad de onda $c = 1\text{ m/s}$, extremos fijos $u(0, t) = u(1, t) = 0$. Perfil inicial en reposo ($g(x) = 0$): pulso triangular con pico en $x = 0.5\text{ m}$ y altura $1.0\text{ m}$:

• $N_x = 4 \implies \Delta x = 0.25\text{ m}$.
• Seleccionamos el número de Courant crítico $C = \frac{c \Delta t}{\Delta x} = 1.0 \implies \Delta t = 0.25\text{ s}$.

Valores en $t = 0$: $u^0 = [0.0, \, 0.5, \, 1.0, \, 0.5, \, 0.0]$

Para el primer paso, al partir del reposo ($g = 0$): $u_i^1 = u_i^0 + \frac{C^2}{2}(u_{i+1}^0 - 2u_i^0 + u_{i-1}^0) = \frac{1}{2}(u_{i+1}^0 + u_{i-1}^0)$

• $u_1^1 = \frac{1}{2}(1.0 + 0.0) = 0.5$
• $u_2^1 = \frac{1}{2}(0.5 + 0.5) = 0.5$
• $u_3^1 = \frac{1}{2}(0.0 + 1.0) = 0.5$ Vector en $t = 0.25\text{ s}$: $u^1 = [0.0, \, 0.5, \, 0.5, \, 0.5, \, 0.0]$.

Para el segundo paso ($t = 0.50\text{ s}$), usamos el esquema completo: $u_i^2 = 2u_i^1 - u_i^0 + (u_{i+1}^1 - 2u_i^1 + u_{i-1}^1) = u_{i+1}^1 + u_{i-1}^1 - u_i^0$

• $u_1^2 = u_2^1 + u_0^1 - u_1^0 = 0.5 + 0.0 - 0.5 = 0.0$
• $u_2^2 = u_3^1 + u_1^1 - u_2^0 = 0.5 + 0.5 - 1.0 = 0.0$
• $u_3^2 = u_4^1 + u_2^1 - u_3^0 = 0.0 + 0.5 - 0.5 = 0.0$ Vector en $t = 0.50\text{ s}$: $u^2 = [0.0, \, 0.0, \, 0.0, \, 0.0, \, 0.0]$. ¡En $t = 0.5\text{ s}$, la onda se ha dividido exactamente en dos paquetes que viajan en direcciones opuestas y se anulan instantáneamente en el centro, reproduciendo de forma exacta el principio de superposición de D'Alembert sin amortiguamiento artificial debido a que $C = 1$!

Errores Comunes en Diferencias Finitas y Resolución de EDPs

1. Violar la Condición de Estabilidad CFL o von Neumann

• Error: Reducir a la mitad el espaciado espacial $\Delta x$ para mejorar la resolución sin ajustar el paso temporal $\Delta t$.
• Consecuencia: En difusión, el parámetro $r = \alpha \Delta t / \Delta x^2$ se cuadruplica. Si supera $0.5$, el factor de amplificación $|G| > 1$ desencadena una explosión numérica inmediata con valores infinitos o NaN.
• Corrección: Recalcular dinámicamente $\Delta t$ cada vez que se modifique la malla espacial para satisfacer $r \le 0.5$ o $C \le 1.0$.

2. Implementar Condiciones de Frontera Neumann con Primer Orden

• Error: Discretizar $\partial u / \partial x = 0$ en el borde como $\frac{u_1 - u_0}{\Delta x} = 0 \implies u_0 = u_1$.
• Consecuencia: Se introduce un error local de orden $\mathcal{O}(\Delta x)$ en la frontera que contamina todo el dominio interior, reduciendo la precisión global del esquema de segundo orden a primer orden.
• Corrección: Emplear puntos fantasma con diferencias centradas de segundo orden $\mathcal{O}(\Delta x^2)$: $u_{-1} = u_1$.

3. Confundir Difusión Física con Difusión Numérica Artificial

• Error: Utilizar esquemas asimétricos hacia adelante o hacia atrás (upwind de primer orden) para advección o transporte de ondas.
• Consecuencia: El término de error de truncamiento actúa como una viscosidad artificial proporcional a $\frac{\Delta x}{2}(1 - C)$, amortiguando picos y ensanchando frentes de onda de forma espuria.
• Corrección: Usar esquemas centrados de segundo orden (como Lax-Wendroff) o esquemas TVD (Total Variation Diminishing).

4. Asumir que Estabilidad Garantiza Convergencia sin Consistencia

• Error: Creer que un algoritmo numérico que no explota numéricamente converge automáticamente a la solución física correcta.
• Consecuencia: De acuerdo con el Teorema de Equivalencia de Lax, la estabilidad garantiza la convergencia si y solo si el esquema es consistente (es decir, el error de truncamiento local tiende a cero cuando $\Delta x, \Delta t \to 0$). Un esquema inestable o inconsistente converge a una solución errónea.
• Corrección: Verificar analíticamente la consistencia mediante series de Taylor antes de implementar cualquier discretización.

5. Ignorar el Manejo de Discontinuidades y Frentes de Choque

• Error: Aplicar diferencias finitas centradas lineales en presencia de ondas de choque o gradientes infinitos.
• Consecuencia: Fenómeno de Gibbs numérico con oscilaciones parásitas severas y sobrepicos espurios antes y después del frente de choque.
• Corrección: Emplear limitadores de flujo o métodos de captura de choque en formulación conservativa.

Aplicaciones Físicas de Diferencias Finitas

1. Conducción Térmica Transitoria en Sólidos y Criogenia

En ingeniería térmica y física del estado sólido, la distribución de calor en piezas estructurales, reactores nucleares o componentes electrónicos se resuelve mediante diferencias finitas 2D/3D: $\rho c_p \frac{\partial T}{\partial t} = \vec{\nabla} \cdot (k \vec{\nabla} T) + Q(\vec{x}, t)$ La discretización permite modelar gradientes térmicos en disipadores de procesadores, blindajes térmicos de reentrada atmosférica y contenedores criogénicos de helio líquido para imanes superconductores.

2. Acústica de Salas y Propagación de Ondas Sísmicas

En geofísica y sismología computacional, la propagación de ondas elásticas $P$ (compresionales) y $S$ (de cizalla) generadas por fallas tectónicas se resuelve con mallas de diferencias finitas en el dominio del tiempo (FDTD): $\rho \frac{\partial^2 \vec{u}}{\partial t^2} = (\lambda + \mu) \vec{\nabla}(\vec{\nabla} \cdot \vec{u}) + \mu \nabla^2 \vec{u}$ Esto permite predecir la aceleración sísmica del terreno y diseñar estructuras sismorresistentes, así como optimizar la respuesta acústica en salas de conciertos.

3. Electromagnetismo Computacional (Método FDTD de Yee)

El algoritmo FDTD formulado por Kane Yee discretiza las ecuaciones de Maxwell en el tiempo entrelazando las componentes espaciales del campo eléctrico $\vec{E}$ y magnético $\vec{H}$: $\vec{\nabla} \times \vec{E} = -\frac{\partial \vec{B}}{\partial t}, \quad \vec{\nabla} \times \vec{H} = \vec{J} + \frac{\partial \vec{D}}{\partial t}$ Es la herramienta reina para diseñar antenas de telefonía móvil 5G, guías de ondas ópticas, metamateriales y sensores nanofotónicos.

4. Ecuación de Schrödinger Dependiente del Tiempo en Mecánica Cuántica

En física cuántica y nanotecnología, la dispersión de paquetes de onda electrónicos a través de barreras de potencial (efecto túnel) se calcula resolviendo: $i \hbar \frac{\partial \psi}{\partial t} = -\frac{\hbar^2}{2m} \frac{\partial^2 \psi}{\partial x^2} + V(x) \psi$ El esquema de Crank-Nicolson se utiliza extensivamente porque garantiza de forma exacta la unitariedad y la conservación de la probabilidad total $\int |\psi|^2 dx = 1$.

5. Dispersión y Transporte de Contaminantes en Hidrodinámica

En física ambiental y oceanografía, el transporte de sustancias disueltas en ríos, acuíferos o atmósfera combina advección hidrodinámica y difusión molecular/turbulenta: $\frac{\partial C}{\partial t} + \vec{v} \cdot \vec{\nabla} C = D \nabla^2 C + S$ La resolución por diferencias finitas permite predecir la pluma de dispersión de fugas químicas accidentales, mareas negras o dispersión de aerosoles contaminantes en ciudades.

Problemas Prácticos: Diferencias Finitas y Estabilidad de Mallas

Problema 1: Inestabilidad Numérica por Violación del Límite de Difusión

Una simulación de difusión de calor 1D con $\alpha = 2.0 \times 10^{-4}\text{ m}^2/\text{s}$ utiliza una malla con espaciado $\Delta x = 0.01\text{ m}$.

1. Calcule el paso de tiempo crítico $\Delta t_{\text{crítico}}$ más allá del cual el esquema FTCS se vuelve inestable.
2. Si un usuario desprevenido ejecuta el código con $\Delta t = 0.3\text{ s}$, demuestre mediante el factor de amplificación de von Neumann qué ocurre con las componentes de frecuencia espacial máxima ($k \Delta x = \pi$).

Solución Paso a Paso:

1. La condición de estabilidad exige: $r = \frac{\alpha \Delta t}{\Delta x^2} \le 0.5$ $\Delta t_{\text{crítico}} = \frac{\Delta x^2}{2 \alpha} = \frac{(0.01)^2}{2(2.0 \times 10^{-4})} = \frac{10^{-4}}{4.0 \times 10^{-4}} = 0.25\text{ s}$
2. Con $\Delta t = 0.3\text{ s}$: $r = \frac{(2.0 \times 10^{-4})(0.3)}{10^{-4}} = 0.60 > 0.50$ El factor de amplificación para el modo Nyquist de la malla ($k \Delta x = \pi$) es: $G = 1 - 4r \sin^2(\pi / 2) = 1 - 4(0.60)(1) = 1 - 2.40 = -1.40$ Como $|G| = 1.40 > 1$, la amplitud de este modo se multiplica por $1.4$ en cada paso temporal, alternando de signo:
   • Tras 10 pasos: $(1.40)^{10} \approx 28.9$ veces mayor.
   • Tras 50 pasos: $(1.40)^{50} \approx 1.8 \times 10^7$ veces mayor. Se genera un patrón espurio en diente de sierra que destruye la simulación en pocas iteraciones.

Problema 2: Implementación de Extremo Aislante con Punto Fantasma

En una barra de difusión térmica de longitud $L = 0.5\text{ m}$, el extremo izquierdo $x = 0$ está térmicamente aislado del ambiente, de modo que el flujo de calor es cero: $\left.\frac{\partial u}{\partial x}\right|_{x=0} = 0$

1. Utilizando una aproximación centrada de segundo orden con punto fantasma $x_{-1} = -\Delta x$, derive la ecuación explícita para la temperatura $u_0^{n+1}$ en el contorno.
2. Si $\Delta x = 0.1\text{ m}$, $\alpha = 1\text{ m}^2/\text{s}$, $\Delta t = 0.004\text{ s}$ y los valores en $t = n$ son $u_0^n = 100^\circ\text{C}$ y $u_1^n = 90^\circ\text{C}$, calcule la temperatura en el extremo en el paso siguiente $u_0^{n+1}$.

Solución Paso a Paso:

1. Condición de Neumann centrada: $\frac{u_1^n - u_{-1}^n}{2 \Delta x} = 0 \implies u_{-1}^n = u_1^n$ Ecuación FTCS para el nodo $i = 0$: $u_0^{n+1} = u_0^n + r (u_1^n - 2u_0^n + u_{-1}^n)$ Sustituyendo $u_{-1}^n = u_1^n$: $u_0^{n+1} = u_0^n + 2r (u_1^n - u_0^n)$
2. Calculamos el parámetro $r$: $r = \frac{\alpha \Delta t}{\Delta x^2} = \frac{(1)(0.004)}{(0.1)^2} = \frac{0.004}{0.01} = 0.40$ Evaluamos: $u_0^{n+1} = 100 + 2(0.40)(90 - 100) = 100 + 0.80(-10) = 100 - 8.0 = 92.0^\circ\text{C}$ La temperatura desciende debido a que el calor fluye hacia el nodo interior adyacente más frío.

Problema 3: Dispersión Numérica en Ondas Acústicas

En la simulación de la ecuación de onda con $c = 340\text{ m/s}$, se usa una celda $\Delta x = 0.1\text{ m}$ y $\Delta t = 0.0002\text{ s}$.

1. Verifique si se satisface la condición de Courant-Friedrichs-Lewy (CFL).
2. Explique por qué trabajar con un número de Courant significativamente menor a 1 (por ejemplo $C = 0.3$) introduce dispersión numérica artificial que deforma los frentes de onda.

Solución Paso a Paso:

1. Número de Courant: $C = \frac{c \Delta t}{\Delta x} = \frac{(340)(0.0002)}{0.1} = \frac{0.068}{0.1} = 0.68$ Dado que $C = 0.68 \le 1.0$, el esquema es formalmente estable.
2. Cuando $C < 1$, la relación de dispersión discreta $\omega_{\text{num}}(k)$ se desvía de la relación física lineal $\omega = c k$: $\sin\left(\frac{\omega_{\text{num}} \Delta t}{2}\right) = C \sin\left(\frac{k \Delta x}{2}\right)$ Como resultado, la velocidad de fase numérica $v_{\text{fase}} = \omega_{\text{num}} / k$ depende del número de onda $k$. Las ondas de longitud de onda corta (alta frecuencia) viajan más lentamente que las de onda larga, provocando oscilaciones espurias (ripples) detrás de frentes de onda originalmente abruptos. Solo cuando $C = 1.0$ se anula idénticamente el error de dispersión en 1D.