Aproximación numérica y error

Fuentes de error numérico: truncamiento, redondeo aritmético de coma flotante, convergencia y estabilidad.

Qué es

La aproximación numérica es la disciplina físico-matemática que modela, cuantifica y controla las discrepancias inevitables que surgen cuando una ley física continua formulada en el continuo real $\mathbb{R}$ se traduce a un algoritmo computable ejecutado en una máquina de precisión finita.

El error total en cualquier cálculo computacional se descompone en tres fuentes fundamentales irreducibles:

1. Error de truncamiento ($E_{\text{trunc}}$): Surge al aproximar un operador analítico continuo infinito (una serie infinita de Taylor, una derivada como límite $\lim_{h \to 0}$, o una integral de Riemann) mediante una fórmula algebraica discreta finita con paso de malla $h$. Su orden asintótico se expresa en notación de Landau $\mathcal{O}(h^p)$, donde $p$ es el orden de convergencia.
2. Error de redondeo ($E_{\text{red}}$): Consecuencia de representar números reales en precisión binaria fija bajo el estándar IEEE-754 (precisión simple de 32 bits con mantisa de 24 bits y $\epsilon_{\text{máq}} \approx 1.19 \times 10^{-7}$; o precisión doble de 64 bits con mantisa de 53 bits y $\epsilon_{\text{máq}} \approx 2.22 \times 10^{-16}$). Cada operación elemental introduce un error relativo acotado por el épsilon de la máquina $\epsilon_{\text{máq}}$.
3. Cancelación catastrófica: Fenómeno patológico que ocurre al restar dos números de magnitud casi idéntica en punto flotante: los dígitos significativos más altos se cancelan entre sí, dejando en los dígitos restantes únicamente el ruido de redondeo previo y perdiendo drásticamente la precisión.

Por qué aparece en física

En física teórica moderna y computacional (astrofísica, mecánica cuántica de muchos cuerpos, dinámica de fluidos geofísicos, física de plasmas), casi ninguna ecuación diferencial no lineal admite soluciones analíticas cerradas:

• Modelar la trayectoria de $N$ cuerpos gravitatorios, el clima terrestre o el confinamiento magnético en un reactor tokamak exige discretizar el espacio y el tiempo.
• Existe una tensión competitiva irreconciliable conocida como el compromiso paso-error (step-size trade-off): al reducir el paso de integración $h$, el error de truncamiento decrece monótonamente como $C_1 h^p$, pero el número de operaciones aritméticas se dispara como $1/h$, acumulando un error de redondeo estocástico que crece como $C_2 \epsilon_{\text{máq}} / h^q$. Por debajo de un paso óptimo $h^*$, disminuir $h$ empeora la precisión global.
• Ignorar la estabilidad numérica y el condicionamiento de matrices físicas conduce a simulaciones físicamente absurdas donde la energía diverge a infinito, se violan leyes de conservación o aparecen inestabilidades numéricas explosivas puramente artificiales.

Idea visual o geométrica

Visualmente, el error total en función del tamaño del paso $h$ se representa en un gráfico logarítmico-logarítmico ($\log E$ frente a $\log h$):

• Para $h$ grande, la curva es una recta con pendiente positiva $+p$ dominada por el truncamiento de Taylor: $\log E \approx p \log h + \text{cte}$.
• Para $h$ minúsculo, la curva quiebra hacia arriba formando una rama con pendiente negativa dominada por la acumulación de redondeo y cancelación numérica.
• El punto más bajo de la curva en forma de "V" define el paso óptimo $h^*$, donde la derivada del error total respecto a $h$ se anula. Calcular con un paso menor a $h^*$ no solo desperdicia tiempo de CPU, sino que produce un resultado numéricamente más erróneo.

Lenguaje simbólico mínimo

1. Error absoluto y relativo: $E_{\text{abs}} = |x_{\text{verdadero}} - x_{\text{aprox}}|, \qquad E_{\text{rel}} = \frac{|x_{\text{verdadero}} - x_{\text{aprox}}|}{|x_{\text{verdadero}}|}$
2. Modelo de punto flotante IEEE-754: $\operatorname{fl}(x) = x (1 + \delta), \qquad |\delta| \leq \epsilon_{\text{máq}}$
3. Cancelación catastrófica (resta de números cercanos $x \approx y$): $\operatorname{fl}(x - y) = (x - y) + x \delta_1 - y \delta_2 = (x - y) \left[ 1 + \frac{x \delta_1 - y \delta_2}{x - y} \right]$ Como $(x - y) \to 0$, el error relativo formal diverge si no se reformula analíticamente la expresión.
4. Número de condición de un problema físico $f(x)$: $\kappa(x) = \left| \frac{x f'(x)}{f(x)} \right|$ Si $\kappa(x) \gg 1$, el problema está mal condicionado (sensible a perturbaciones infinitesimales).

Qué significa el resultado

Comprender el error numérico transforma un código de simulación física en una herramienta científica rigurosa:

• Un resultado numérico carece de valor si no va acompañado de una cota de error comprobable o una demostración de convergencia asintótica en malla.
• En simulaciones orbitales o mecánicas, verificar que el error global decrece a la tasa teórica predicha ($h^2, h^4$) garantiza que no hay errores lógicos en el código ni inestabilidades numéricas ocultas.

Cuándo no debes usarlo

1. No disminuyas ciegamente el paso $h$ sin verificar el límite de redondeo: En derivadas numéricas simples por diferencias finitas en precisión doble, si eliges $h = 10^{-16}$, el resultado no tenderá a la derivada real, sino que arrojará exactamente cero o números aleatorios debido a la cancelación catastrófica.
2. No utilices precisión simple (Float32) para integraciones dinámicas de largo plazo: La acumulación sistemática de redondeo ($10^{-7}$) tras millones de pasos de integración destruye la conservación de la energía en simulaciones moleculares o gravitatorias.
3. No confundas aproximación numérica con simplificación de la física: El error numérico mide la discrepancia entre el algoritmo computacional y la ecuación diferencial matemática elegida; no compensa ni detecta si el modelo físico de partida omitió rozamientos o efectos relativistas fundamentales.

Resumen operativo

• Identifica la formulación matemática continua y el algoritmo de discretización.
• Determina el orden teórico del error de truncamiento $\mathcal{O}(h^p)$.
• Analiza si la fórmula contiene restas de magnitudes casi idénticas o raíces cuadradas cercanas que induzcan cancelación catastrófica; si es así, reescribe algebraicamente multiplicando por el conjugado o usando series de Taylor analíticas.
• Selecciona el paso de malla $h$ equilibrando el truncamiento frente al límite de precisión del formato IEEE-754.
• Realiza siempre un ensayo de convergencia en malla (grid convergence test) reduciendo $h$ a la mitad ($h, h/2, h/4$) y calculando el cociente de convergencia experimental.

Procedimiento de Control y Mitigación de Errores Numéricos

Antes de calcular

Antes de ejecutar una rutina computacional en física:

• Identifica la precisión de punto flotante del entorno de cómputo (Float32: $\epsilon_{\text{máq}} \approx 1.2 \times 10^{-7}$, 7 dígitos decimales; Float64: $\epsilon_{\text{máq}} \approx 2.2 \times 10^{-16}$, 16 dígitos decimales).
• Revisa las expresiones algebraicas en busca de operaciones vulnerables a cancelación catastrófica (por ejemplo, $1 - \cos(x)$ para $x \to 0$, o $\sqrt{x + 1} - \sqrt{x}$ para $x \gg 1$).

Paso 1: identifica el objeto matemático

Determina la naturaleza del algoritmo numérico:

• Derivada numérica por diferencias finitas: Truncamiento $\mathcal{O}(h)$ en diferencias progresivas o $\mathcal{O}(h^2)$ en diferencias centradas.
• Integración numérica (cuadratura): Truncamiento $\mathcal{O}(h^2)$ en regla trapezoidal o $\mathcal{O}(h^4)$ en regla de Simpson.
• Resolución de EDOs (Runge-Kutta): Error de truncamiento local $\mathcal{O}(h^{p+1})$ y global $\mathcal{O}(h^p)$.

Paso 2: elige la operación (Reformulación contra cancelación)

Si detectas resta de números casi idénticos, reescribe la fórmula analíticamente antes de programarla:

• Resta de raíces: $\sqrt{a} - \sqrt{b} = \frac{a - b}{\sqrt{a} + \sqrt{b}}$ (multiplicar por el conjugado elimina la cancelación).
• Trigonometría para ángulos pequeños: $1 - \cos(x) = 2 \sin^2(x/2)$ o su expansión de Taylor $\frac{x^2}{2} - \frac{x^4}{24} + \dots$.
• Ecuaciones cuadráticas con $b^2 \gg 4 a c$: $x_1 = \frac{-b - \operatorname{sgn}(b)\sqrt{b^2 - 4ac}}{2a}, \qquad x_2 = \frac{2c}{-b - \operatorname{sgn}(b)\sqrt{b^2 - 4ac}} = \frac{c}{a x_1}$

Paso 3: conserva unidades y signos

• Los errores absolutos $E_{\text{abs}}$ mantienen estrictamente las unidades de la magnitud física calculada.
• Los errores relativos $E_{\text{rel}}$ y el épsilon de la máquina $\epsilon_{\text{máq}}$ son números puros adimensionales.
• Al estimar derivadas numéricas $f'(x) \approx \frac{f(x+h) - f(x-h)}{2h}$, el paso $h$ debe tener las mismas unidades físicas que la variable independiente $x$.

Paso 4: calcula (Determinación del paso óptimo $h^*$)

Para un método con error de truncamiento $\alpha h^p$ y error de redondeo $\frac{\epsilon_{\text{máq}}}{h^q}$:

1. Plantea el error total: $E_{\text{total}}(h) = \alpha h^p + \beta \frac{\epsilon_{\text{máq}}}{h^q}$
2. Deriva respecto a $h$ e iguala a cero: $\frac{d E_{\text{total}}}{dh} = p \alpha h^{p-1} - q \beta \frac{\epsilon_{\text{máq}}}{h^{q+1}} = 0$
3. Despeja el paso óptimo: $h^* = \left( \frac{q \beta \epsilon_{\text{máq}}}{p \alpha} \right)^{\frac{1}{p + q}}$ (Por ejemplo, para diferencias centradas con $p=2, q=1$: $h^* \sim \epsilon_{\text{máq}}^{1/3}$. En Float64, $h^* \sim (10^{-16})^{1/3} \approx 10^{-5} - 10^{-6}$).

Paso 5: interpreta el resultado (Test de convergencia de Richardson)

Ejecuta el cálculo con tres pasos sucesivamente divididos por dos ($h, h/2, h/4$) obteniendo las soluciones aproximadas $A_1, A_2, A_3$:

1. Calcula el cociente de convergencia: $R = \frac{A_1 - A_2}{A_2 - A_3} \approx 2^p$
2. Si $R \approx 2^p$ (por ejemplo, $R \approx 4$ para métodos de segundo orden $p=2$, o $R \approx 16$ para cuarto orden $p=4$), el algoritmo está convergiendo limpiamente en el régimen asintótico sin contaminación de redondeo.
3. Extrapolación de Richardson para eliminar el error dominante: $A_{\text{mejorado}} = A_3 + \frac{A_3 - A_2}{2^p - 1}$

Checklist de autocontrol

• ¿He verificado que el paso $h$ elegido sea mayor que el paso óptimo de redondeo $h^*$?
• ¿He eliminado todas las cancelaciones catastróficas mediante conjugados o expansiones analíticas?
• ¿Es el cociente de convergencia experimental $R$ coherente con el orden teórico $2^p$?
• ¿He utilizado Float64 para cálculos que involucran miles o millones de iteraciones temporales?
• ¿Se conservan invariantes físicos fundamentales (energía, momento lineal) dentro de la tolerancia de truncamiento?

Ejemplo guiado: Cancelación catastrófica en cinemática relativista

Situación física

Una partícula elemental de masa en reposo $m = 1.000 \text{ GeV}/c^2$ se acelera en el Gran Colisionador de Hadrones (LHC) hasta alcanzar una energía total ultra-relativista $E = 1000.000000000100 \text{ GeV}$. Se desea calcular computacionalmente el momento lineal relativista $p$ utilizando la relación de Einstein: $E^2 = p^2 c^2 + m^2 c^4 \implies p c = \sqrt{E^2 - m^2 c^4}$ En unidades naturales ($c=1$), $p = \sqrt{E^2 - m^2}$. Como $E \gg m$, el momento $p$ es prácticamente igual a $E$. Si un programador calcula directamente la energía cinética relativista como: $E_k = E - m$ y la diferencia de velocidades respecto a la luz $\Delta v = c - v = 1 - \frac{p}{E} = 1 - \sqrt{1 - \frac{m^2}{E^2}}$, se produce una pérdida catastrófica de dígitos significativos en precisión simple e incluso en precisión doble estándar.

Traducción matemática

Evaluamos la expresión para $E = 10^4 \text{ GeV}$ y $m = 1 \text{ GeV}$: $\gamma = \frac{E}{m} = 10^4, \qquad \frac{m^2}{E^2} = 10^{-8}$ La fórmula directa calcula: $1 - \sqrt{1 - 10^{-8}}$ En punto flotante con precisión limitada, $\sqrt{1 - 10^{-8}} = 0.9999999950000000\dots$ Al restar $1 - 0.999999995\dots$, los primeros 8 dígitos idénticos se cancelan, perdiendo la mitad de la mantisa en Float64 y el $100\%$ de los dígitos en Float32 (donde $1 - 10^{-8}$ redondea directamente a $1.0$, dando $\Delta v = 0$).

Cálculo paso a paso

1. Reformulación analítica mediante multiplicación por el conjugado: $\Delta v = 1 - \sqrt{1 - \frac{m^2}{E^2}} = \frac{\left(1 - \sqrt{1 - \frac{m^2}{E^2}}\right)\left(1 + \sqrt{1 - \frac{m^2}{E^2}}\right)}{1 + \sqrt{1 - \frac{m^2}{E^2}}} = \frac{1 - \left(1 - \frac{m^2}{E^2}\right)}{1 + \sqrt{1 - \frac{m^2}{E^2}}} = \frac{m^2 / E^2}{1 + \sqrt{1 - \frac{m^2}{E^2}}}$
2. Evaluación numérica sin cancelación: El denominador ahora es una suma de números positivos: $1 + \sqrt{1 - 10^{-8}} \approx 1 + 0.999999995 = 1.999999995$ El numerador es: $\frac{m^2}{E^2} = 10^{-8}$ Por lo tanto: $\Delta v = \frac{10^{-8}}{1.999999995} \approx 5.0000000125 \times 10^{-9} \text{ (en unidades de } c)$
3. Alternativa por expansión binomial de Taylor: $\sqrt{1 - \epsilon} = 1 - \frac{1}{2}\epsilon - \frac{1}{8}\epsilon^2 - \dots \quad \text{con } \epsilon = \frac{m^2}{E^2} = 10^{-8}$ $1 - \sqrt{1 - \epsilon} = \frac{1}{2}\epsilon + \frac{1}{8}\epsilon^2 = \frac{1}{2}(10^{-8}) + \frac{1}{8}(10^{-16}) = 5.0000000000000125 \times 10^{-9}$

Lectura del resultado

La reformulación algebraica preserva la totalidad de los 16 dígitos significativos de Float64. La partícula se desplaza a una velocidad inferior a la de la luz por apenas $\Delta v \approx 5.0 \times 10^{-9} c \approx 1.5 \text{ m/s}$. El cálculo directo ingenuo habría devuelto 0 m/s o un valor con un error relativo del $50\%$.

Error común evitado

Calcular diferencias de raíces cuadradas directamente en software sin comprobar si el argumento relativo es minúsculo. El uso de identidades de conjugados elimina de raíz la cancelación catastrófica sin coste computacional adicional.

Ejemplo inverso: Búsqueda del paso óptimo en diferenciación numérica

Situación física

Se desea calcular numéricamente la aceleración instantánea $a(t) = v'(t)$ de un cohete espacial a partir de mediciones de velocidad $v(t) = \sin(t)$ en $t = 1.0 \text{ s}$ utilizando el método de diferencias centradas en una máquina con precisión estándar Float64 ($\epsilon_{\text{máq}} \approx 2.22 \times 10^{-16}$). La derivada teórica exacta es $v'(1.0) = \cos(1.0) \approx 0.5403023058681398$.

Resultado dado

Se realizan pruebas numéricas con pasos decrecientes $h$ obteniendo los siguientes errores absolutos:

• $h = 10^{-2}$: $E_{\text{abs}} = 9.00 \times 10^{-6}$ (dominado por truncamiento $\propto h^2$).
• $h = 10^{-4}$: $E_{\text{abs}} = 9.00 \times 10^{-10}$ (truncamiento).
• $h = 10^{-6}$: $E_{\text{abs}} = 2.11 \times 10^{-11}$ (punto mínimo del error).
• $h = 10^{-9}$: $E_{\text{abs}} = 8.35 \times 10^{-8}$ (el error empeora drásticamente).
• $h = 10^{-14}$: $E_{\text{abs}} = 5.40 \times 10^{-2}$ (pérdida total de precisión).

Reconstrucción analítica

La fórmula de diferencias centradas es: $f'(x) \approx \frac{f(x+h) - f(x-h)}{2h}$ El error de truncamiento por serie de Taylor es: $E_{\text{trunc}} = \frac{h^2}{6} |f'''(x)| = \frac{h^2}{6} |-\cos(1.0)| \approx 0.09005 h^2$ El error de redondeo al evaluar el numerador en punto flotante es aproximadamente $\epsilon_{\text{máq}} |f(x)|$, el cual al dividirse por $2h$ genera: $E_{\text{red}} \approx \frac{\epsilon_{\text{máq}}}{h} |\sin(1.0)| \approx 0.8415 \frac{\epsilon_{\text{máq}}}{h}$ Minimizando el error total $E(h) = 0.09005 h^2 + 0.8415 \frac{\epsilon_{\text{máq}}}{h}$: $\frac{dE}{dh} = 2(0.09005) h - 0.8415 \frac{\epsilon_{\text{máq}}}{h^2} = 0 \implies h^{*3} = \frac{0.8415 \epsilon_{\text{máq}}}{0.1801} \approx 4.67 \epsilon_{\text{máq}}$ Para $\epsilon_{\text{máq}} \approx 2.22 \times 10^{-16}$: $h^* = (4.67 \times 2.22 \times 10^{-16})^{1/3} = (1.037 \times 10^{-15})^{1/3} \approx 1.01 \times 10^{-5}$

Qué pregunta física permite responder

¿Por qué elegir $h = 10^{-14}$ en vez de $h = 10^{-5}$ destruye el cálculo de la aceleración del cohete? Porque en $h = 10^{-14}$, los valores $v(t+h)$ y $v(t-h)$ coinciden en casi todos sus bits binarios. Su resta produce cancelación catastrófica, dejando un residuo de ruido puro en el último bit que se multiplica por $1/(2h) = 10^{14} / 2$, amplificando el ruido de punto flotante por catorce órdenes de magnitud.

Verificación de límites

Cuando $h \to 0$, $E_{\text{trunc}} \to 0$ pero $E_{\text{red}} \to \infty$. Solo en una máquina hipotética con infinitos bits de mantisa ($\epsilon_{\text{máq}} = 0$) sería óptimo tomar el límite infinitesimal $h \to 0$.

Errores comunes en Aproximación Numérica y Error

1. Disminuir el paso de integración $h$ creyendo que el error siempre disminuye

Parece razonable porque la teoría de cálculo infinitesimal define las derivadas e integrales en el límite $h \to 0$. Falla porque en una computadora real con precisión finita (IEEE-754), el error total tiene forma de "V": al disminuir $h$, el error de truncamiento decrece pero el error de redondeo crece como $\epsilon_{\text{máq}} / h^q$. Por debajo del paso óptimo $h^*$, disminuir $h$ multiplica el error global. Se detecta cuando al reducir $h$ el error experimental comienza a oscilar o crecer abruptamente en lugar de decaer como $h^p$. Se corrige calculando el paso óptimo $h^* \sim \epsilon_{\text{máq}}^{1/(p+q)}$ y realizando análisis de convergencia en malla. Mini-ejemplo: calcular $f'(x) \approx \frac{f(x+h) - f(x)}{h}$ con $h = 10^{-16}$ en Float64 devuelve $0.0$ en vez de la derivada real, porque en el numerador $f(x+10^{-16})$ redondea a $f(x)$.

2. Restar números de magnitud casi idéntica en punto flotante (cancelación catastrófica)

Parece razonable porque las computadoras son máquinas deterministas de alta precisión con 16 dígitos decimales. Falla porque si dos números coinciden en sus primeros 14 dígitos decimales, su resta elimina esos 14 dígitos significativos, dejando solo 2 dígitos de información real seguidos de ruido de redondeo en los bits inferiores. Se detecta si aparecen resultados numéricos con pérdidas súbitas de precisión en diferencias de raíces cuadradas, funciones trigonométricas en ángulos pequeños o polinomios cuadráticos. Se corrige reescribiendo la expresión analíticamente mediante identidades algebraicas (multiplicar por conjugados, expansiones de Taylor o identidades de ángulo mitad). Mini-ejemplo: en precisión Float32, evaluar $1 - \cos(10^{-4})$ da $0.0$ directamente; usar $2 \sin^2(10^{-4}/2)$ devuelve con total exactitud $5.000000 \times 10^{-9}$.

3. Comparar números de punto flotante usando igualdad exacta

Parece razonable porque matemáticamente dos expresiones analíticamente idénticas deberían ser iguales. Falla porque las operaciones en punto flotante no son asociativas ni distributivas en binario: $(0.1 + 0.2) + 0.3 \neq 0.1 + (0.2 + 0.3)$ en IEEE-754. La expresión 0.1 + 0.2 == 0.3 evalúa a falso en la inmensa mayoría de lenguajes de programación. Se detecta si un bucle iterativo o criterio de parada no termina nunca o se cuelga en un bucle infinito. Se corrige comparando siempre mediante tolerancias absolutas y relativas: abs(x - y) <= atol + rtol * abs(y). Mini-ejemplo: un algoritmo que itera mientras x != 1.0 con x += 0.1 entrará en bucle infinito porque x nunca valdrá exactamente 1.0000000000000000.

4. Confundir estabilidad numérica con precisión del método

Parece razonable porque un método de orden superior (como Runge-Kutta 4) es más preciso que uno de orden bajo (como Euler). Falla porque un método puede ser formalmente muy preciso pero inestable si el problema físico es "rígido" (stiff) o si el paso temporal supera el radio de estabilidad del algoritmo. Se detecta si la solución numérica explota bruscamente hacia $\pm \infty$ tras haber parecido converger durante los primeros pasos temporales. Se corrige analizando la región de estabilidad lineal absoluta en el plano complejo $z = \lambda \Delta t$ y adoptando métodos implícitos (como Crank-Nicolson o Radau) para ecuaciones diferenciales rígidas. Mini-ejemplo: en cinética química con constantes de reacción dispares ($k_1 = 1, k_2 = 10^6$), usar RK4 explícito exige un $\Delta t < 10^{-6}$ no por precisión, sino para evitar que la solución explote numéricamente.

5. Acumular sumas de series en orden descendente de magnitud

Parece razonable porque sumar términos grandes primero parece capturar "lo más importante" al inicio. Falla porque al sumar un término minúsculo a un acumulador gigante, los bits del término pequeño se desplazan fuera de la mantisa de punto flotante y se pierden por completo (absorción numérica). Se detecta si al sumar $N = 10^7$ números positivos pequeños su suma se estanca y no crece a partir de cierto umbral. Se corrige ordenando los términos en orden ascendente (de menor a mayor valor absoluto) o utilizando el algoritmo de compensación de suma de Kahan. Mini-ejemplo: en Float32, si un acumulador vale $10^8$, sumarle $1.0$ no produce ningún efecto: $10^8 + 1.0 = 10^8$. Sumar un millón de veces $1.0$ al acumulador dejará el total en $10^8$ sin cambiar, perdiendo un millón de unidades.

Regla rápida de detección

• Si estás restando términos: comprueba si su cociente relativo está entre $0.99$ y $1.01$. Si es así, hay cancelación catastrófica.
• Si reduces el paso $h$: grafica el error en escala log-log. Si la curva deja de bajar y sube, has cruzado el límite de redondeo.

Aplicaciones físicas de Aproximación Numérica y Error

1. Astrodinámica: cálculo de trayectorias interplanetarias y semiejes orbitales

Aparece en ingeniería aeroespacial al integrar las ecuaciones de movimiento de sondas espaciales interplanetarias (como las misiones Voyager o New Horizons) bajo la influencia gravitacional del Sol y los planetas del Sistema Solar. Se calcula integrando las ecuaciones diferenciales ordinarias $\ddot{\mathbf{r}} = -\sum G M_i \frac{\mathbf{r} - \mathbf{r}_i}{|\mathbf{r} - \mathbf{r}_i|^3}$ durante periodos de décadas temporales. El resultado da la posición y velocidad de la sonda para maniobras de asistencia gravitacional (gravity assist) con márgenes de error inferiores a unas pocas decenas de kilómetros en encuentros a miles de millones de kilómetros de la Tierra. El error peligroso es la acumulación de error de redondeo y la disipación artificial de energía inducida por métodos no conservativos (como RK4 estándar); en una integración de 30 años con pasos de 1 hora, un error de redondeo relativo de $10^{-15}$ por paso se acumula sistemáticamente deformando el semieje mayor orbital. Como control final se utilizan integradores simplécticos de orden superior (Gauss-Radau o Yoshida) y se monitorea la invariancia estricta del vector de Runge-Lenz y la energía mecánica total.

2. Mecánica cuántica computacional: resolución de la ecuación de Schrödinger estacionaria

Aparece en física atómica y del estado sólido al calcular los autoestados de energía y funciones de onda $\psi(x)$ de pozos de potencial cuánticos no armónicos o moléculas diatómicas mediante el método de disparo (shooting method) o discretización por diferencias finitas de la ecuación $-\frac{\hbar^2}{2m} \frac{d^2\psi}{dx^2} + V(x)\psi = E\psi$. Se calcula aproximando la segunda derivada mediante el algoritmo de Numerov de sexto orden $\mathcal{O}(h^6)$. El resultado da los niveles de energía discretos cuantizados $E_n$ y la probabilidad de tunelamiento cuántico en barreras de potencial de semiconductores. El error peligroso es elegir un paso de malla espacial $h$ demasiado grande en las regiones clásicamente prohibidas donde la función de onda decae exponencialmente; si el truncamiento numérico no resuelve la curvatura exponencial, la solución diverge asintóticamente a $\pm \infty$ impidiendo encontrar la energía propia. Como control final se verifica la condición de ortonormalidad de las funciones de onda $\int |\psi_n|^2 dx = 1$ y se contrastan los autovalores con cálculos de perturbaciones analíticas de primer orden.

3. Dinámica de fluidos geofísicos: modelos de predicción meteorológica y clima

Aparece en meteorología numérica al resolver las ecuaciones primitivas hidrodinámicas de Navier-Stokes en la atmósfera terrestre en mallas globales tridimensionales. Se calcula discretizando las derivadas espaciales de presión y advección de viento mediante esquemas en diferencias finitas y transformadas espectrales esféricas. El resultado da la evolución temporal de frentes térmicos, borrascas y huracanes con horizontes de predicción de varios días. El error peligroso es violar la condición de Courant-Friedrichs-Lewy (CFL: $C = \frac{u \Delta t}{\Delta x} \leq 1$); si el viento $u$ o una onda de gravedad atmosférica viaja más rápido que una celda de malla $\Delta x$ por paso temporal $\Delta t$, el error de truncamiento se amplifica exponencialmente generando una explosión numérica no física en cuestión de pocos pasos de simulación. Como control final se implementan pasos temporales adaptativos gobernados por el número de CFL local y se añade viscosidad numérica controlada (hiperdifusión) para disipar vórtices por debajo de la resolución de la malla sin alterar las ondas planetarias de Rossby.

4. Física del plasma y fusión nuclear: seguimiento de partículas en reactores Tokamak

Aparece en física de plasmas en experimentos de fusión magnética (como ITER) al calcular las órbitas de deriva de iones y partículas alfa confinadas por campos magnéticos helicoidales toroidales intensos. Se calcula integrando la fuerza de Lorentz $\mathbf{F} = q (\mathbf{E} + \mathbf{v} \times \mathbf{B})$ mediante el algoritmo cinemático de Boris. El resultado da el radio de Larmor del movimiento giroscópico y verifica si las partículas confinadas impactan contra la pared del reactor (pérdida de confinamiento). El error peligroso es emplear integradores estándar de Runge-Kutta para el término de rotación magnética $\mathbf{v} \times \mathbf{B}$; los métodos estándar no conservan el momento magnético adiabático $\mu = m v_\perp^2 / (2B)$ ni el radio de giro, haciendo que las partículas ganen energía cinética numérica artificial y escapen del campo magnético en simulaciones largas. Como control final el integrador de Boris preserva exactamente la magnitud de la velocidad en campos magnéticos puros ($|\mathbf{v}_{n+1}| = |\mathbf{v}_n|$), garantizando la conservación absoluta de la energía cinética sin importar la duración de la simulación.

5. Óptica electromagnética: simulación FDTD de metamateriales y guías de ondas

Aparece en fotónica y telecomunicaciones al simular la propagación de pulsos de luz láser a través de cristales fotónicos y guías dieléctricas resolviendo las ecuaciones de Maxwell en el dominio del tiempo mediante el método de diferencias finitas en el dominio del tiempo (FDTD / malla de Yee). Se calcula entrelazando las componentes espaciales y temporales del campo eléctrico $\mathbf{E}$ y magnético $\mathbf{H}$ en una cuadrícula cúbica desplazada en medio paso de tiempo y espacio. El resultado da la transmisión espectral, los diagramas de bandas fotónicas y la dispersión cromática de fibras ópticas de última generación. El error peligroso es la dispersión numérica: en una malla discreta, la velocidad de propagación de las ondas electromagnéticas simuladas $v_{\text{num}}$ no es estrictamente igual a $c$, sino que varía artificialmente en función de la dirección de propagación en la cuadrícula y de la longitud de onda en relación a $\Delta x$. Si la malla es gruesa, un pulso óptico ultracorto se ensancha y deforma artificialmente debido al error de truncamiento de la discretización espacial. Como control final se fija la resolución espacial a un mínimo de 20 celdas por longitud de onda en el dieléctrico ($\Delta x \leq \lambda_{\text{mín}} / (20 n_{\text{ref}})$) y se verifica que el número de Courant se ajuste al límite óptimo de estabilidad multidimensional.

Práctica guiada: Aproximación Numérica y Error

Nivel 1: reconocimiento

Para cada una de las siguientes operaciones numéricas, identifica si el error predominante es de truncamiento, de redondeo de punto flotante, o de cancelación catastrófica:

1. Calcular el valor de $\pi$ sumando los primeros 100 términos de la serie de Leibniz $\sum_{k=0}^{99} \frac{(-1)^k}{2k+1}$.
2. Evaluar $y = \sqrt{1000001} - 1000$ en una calculadora de precisión fija.
3. Evaluar $f(x) = \frac{1}{3} x$ almacenando el resultado como Float32 en memoria.
4. Integrar una órbita planetaria newtoniana durante 10 millones de pasos temporales acumulando una deriva constante en el semieje mayor.

Solución:

1. Error de truncamiento (se corta una serie infinita en un número finito de términos).
2. Cancelación catastrófica (resta de dos números de magnitud casi idéntica $\approx 1000.0005 - 1000$).
3. Error de redondeo (el número $1/3$ es una fracción periódica binaria infinita $0.010101\dots_2$ que debe truncarse a la mantisa de 24 bits).
4. Error de redondeo acumulado estocástico/sistemático (drift numérico por propagación de punto flotante a largo plazo).

Nivel 2: aplicación directa

Se aproxima la primera derivada de $f(x) = e^x$ en $x=0$ mediante la fórmula de diferencias progresivas hacia adelante: $D_+(h) = \frac{e^h - 1}{h}$ El valor analítico exacto es $f'(0) = e^0 = 1$. Calcula analíticamente el error de truncamiento en función de $h$ mediante la serie de Taylor de $e^h$ y evalúa numéricamente el error para $h = 0.1$ y $h = 0.01$.

Solución:

1. Expansión de Taylor: $e^h = 1 + h + \frac{h^2}{2} + \frac{h^3}{6} + \mathcal{O}(h^4)$
2. Sustitución en el operador: $D_+(h) = \frac{(1 + h + \frac{h^2}{2} + \frac{h^3}{6} + \dots) - 1}{h} = \frac{h + \frac{h^2}{2} + \frac{h^3}{6} + \dots}{h} = 1 + \frac{h}{2} + \frac{h^2}{6} + \dots$
3. Error de truncamiento analítico: $E_{\text{trunc}} = |D_+(h) - 1| = \frac{h}{2} + \mathcal{O}(h^2)$ El método es de primer orden $\mathcal{O}(h)$.
4. Evaluación numérica:
   • Para $h = 0.1$: $E_{\text{trunc}} \approx \frac{0.1}{2} = 0.05$ (error del $5\%$; valor numérico $D_+(0.1) = \frac{e^{0.1} - 1}{0.1} \approx \frac{1.105171 - 1}{0.1} = 1.05171$).
   • Para $h = 0.01$: $E_{\text{trunc}} \approx \frac{0.01}{2} = 0.005$ (error del $0.5\%$; valor numérico $D_+(0.01) \approx 1.005017$). Al reducir $h$ en un factor 10, el error se reduce exactamente en un factor 10 (convergencia lineal).

Nivel 3: elección de herramienta

Necesitas calcular en un software de astrofísica la energía potencial gravitatoria diferencial entre dos puntos situados a distancias $r_1$ y $r_2$ muy próximas entre sí ($r_2 = r_1 + \delta$ con $\delta \ll r_1$): $\Delta U = -G M m \left( \frac{1}{r_2} - \frac{1}{r_1} \right)$ ¿Cómo debes implementar algebraicamente esta fórmula para garantizar estabilidad numérica absoluta sin cancelación?

Solución: Evaluar directamente $\frac{1}{r_2} - \frac{1}{r_1}$ produce cancelación catastrófica de dígitos. Se debe operar algebraicamente antes de computar: $\frac{1}{r_2} - \frac{1}{r_1} = \frac{r_1 - r_2}{r_1 r_2} = \frac{-\delta}{r_1 (r_1 + \delta)}$ Sustituyendo en la energía potencial: $\Delta U = -G M m \left( \frac{-\delta}{r_1 (r_1 + \delta)} \right) = \frac{G M m \delta}{r_1 (r_1 + \delta)}$ Esta forma matemática involucra únicamente productos y sumas de números positivos. Elimina por completo la resta de términos cercanos, preservando la totalidad de los 16 dígitos significativos de Float64 para cualquier $\delta > 0$, incluso cuando $\delta = 10^{-12} r_1$.

Nivel 4: interpretación física

Un físico computacional ejecuta una simulación del oscilador armónico simple $m \ddot{x} + k x = 0$ con el método de Euler explícito ($x_{n+1} = x_n + v_n \Delta t$; $v_{n+1} = v_n - \omega^2 x_n \Delta t$). Tras 1000 oscilaciones observa que la amplitud del oscilador ha crecido exponencialmente y la energía mecánica se ha multiplicado por 10. Explica por qué este fenómeno no es un error de redondeo ni una física real, sino un error de estabilidad del método numérico.

Solución: Calculamos la matriz de amplificación del método de Euler explícito: $\begin{pmatrix} x_{n+1} \\ v_{n+1} \end{pmatrix} = \begin{pmatrix} 1 & \Delta t \\ -\omega^2 \Delta t & 1 \end{pmatrix} \begin{pmatrix} x_n \\ v_n \end{pmatrix}$ El determinante de la matriz de transición es: $\det(A) = 1 \cdot 1 - (\Delta t)(-\omega^2 \Delta t) = 1 + \omega^2 \Delta t^2$ El determinante representa el factor de amplificación del volumen en el espacio fásico $(x, v)$ en cada paso temporal. Dado que $\det(A) = 1 + \omega^2 \Delta t^2 > 1$ para cualquier $\Delta t > 0$, los autovalores de la matriz tienen módulo estrictamente mayor que 1: $|\lambda| = \sqrt{1 + \omega^2 \Delta t^2} > 1$. Interpretación física: El método de Euler explícito introduce artificialmente una "fuerza anti-amortiguadora" numérica que bombea energía espuria al sistema en cada paso temporal. No es un error de redondeo; es una inestabilidad algorítmica intrínseca. Para sistemas conservativos deben emplearse métodos simplécticos (Euler simpléctico o Verlet) donde $\det(A) \equiv 1$.

Nivel 5: error intencionado

Para resolver la ecuación cuadrática $a x^2 + b x + c = 0$ con $a = 1$, $b = 10^8$, $c = 1$, un estudiante utiliza la fórmula cuadrática escolar estándar: $x_1 = \frac{-b + \sqrt{b^2 - 4ac}}{2a}, \qquad x_2 = \frac{-b - \sqrt{b^2 - 4ac}}{2a}$ El software reporta $x_1 = 0.0$ y $x_2 = -10^8$. Demuestra por qué $x_1 = 0$ es matemáticamente falso y calcula el valor correcto utilizando fórmulas numéricamente estables.

Fallo conceptual: Como $b = 10^8$, $b^2 = 10^{16}$. Al calcular $b^2 - 4ac = 10^{16} - 4$, en precisión doble Float64 el número 4 está por debajo de los 16 dígitos de la mantisa respecto a $10^{16}$, por lo que $\sqrt{b^2 - 4ac}$ redondea exactamente a $10^8$. Al evaluar $x_1 = \frac{-10^8 + 10^8}{2} = \frac{0}{2} = 0.0$, la resta de los dos números idénticos destruye completamente la raíz más pequeña. Sin embargo, por el teorema de Cardano-Vieta sabemos que el producto de las raíces es $x_1 \cdot x_2 = c/a = 1/1 = 1$. Si $x_2 \approx -10^8$, entonces la raíz real debe ser $x_1 = 1 / x_2 = -10^{-8} \neq 0$.

Corrección numéricamente estable: Calculamos primero la raíz sin cancelación (usando el mismo signo que $-b$): $x_2 = \frac{-b - \operatorname{sgn}(b)\sqrt{b^2 - 4ac}}{2a} = \frac{-10^8 - 10^8}{2} = -1.0 \times 10^8$ Y obtenemos la segunda raíz mediante la relación de Vieta: $x_1 = \frac{c}{a x_2} = \frac{1}{1 \cdot (-10^8)} = -1.0 \times 10^{-8}$ Sustituyendo $x_1 = -10^{-8}$ en la ecuación original: $(-10^{-8})^2 + 10^8(-10^{-8}) + 1 = 10^{-16} - 1 + 1 = 10^{-16} \approx 0$ La raíz correcta es exactamente $-10^{-8}$, mientras que el cálculo estándar arrojó un error relativo del $100\%$ ($0.0$).