Derivación numérica

Aproximación de derivadas mediante esquemas de diferencias finitas progresivas, regresivas y centradas.

Qué es

La derivación numérica es la aproximación computacional de la tasa instantánea de cambio o derivada de una función $f'(x)$ (o derivadas de orden superior como la aceleración $f''(x)$) cuando la función se conoce únicamente en una malla discreta de puntos experimentales o cuando su forma analítica es demasiado intrincada para derivar simbólicamente.

Se fundamenta en el desarrollo en serie de Taylor alrededor del punto de evaluación $x$:

1. Diferencia progresiva hacia adelante (Forward difference): $f'(x) \approx \frac{f(x+h) - f(x)}{h}, \qquad E_{\text{trunc}} = -\frac{h}{2} f''(\xi) = \mathcal{O}(h)$
2. Diferencia regresiva hacia atrás (Backward difference): $f'(x) \approx \frac{f(x) - f(x-h)}{h}, \qquad E_{\text{trunc}} = +\frac{h}{2} f''(\xi) = \mathcal{O}(h)$
3. Diferencia centrada (Central difference): $f'(x) \approx \frac{f(x+h) - f(x-h)}{2h}, \qquad E_{\text{trunc}} = -\frac{h^2}{6} f'''(\xi) = \mathcal{O}(h^2)$ Por simetría par, los errores lineales de primer orden se cancelan idénticamente, elevando la aproximación a segundo orden $\mathcal{O}(h^2)$.
4. Segunda derivada centrada (operador laplaciano 1D): $f''(x) \approx \frac{f(x+h) - 2f(x) + f(x-h)}{h^2}, \qquad E_{\text{trunc}} = -\frac{h^2}{12} f^{(4)}(\xi) = \mathcal{O}(h^2)$

Por qué aparece en física

En física teórica y experimental, las leyes dinámicas fundamentales son ecuaciones diferenciales que conectan magnitudes con sus derivadas temporales o espaciales:

• Cinemática y dinámica newtoniana: La velocidad instantánea es la primera derivada temporal de la posición $v(t) = \dot{x}(t)$ y la aceleración es la segunda derivada $a(t) = \ddot{x}(t) = [x(t+\Delta t) - 2x(t) + x(t-\Delta t)] / \Delta t^2$.
• Mecánica analítica y campos conservativos: Las fuerzas físicas conservativas se obtienen como el gradiente negativo de la energía potencial: $\mathbf{F}(\mathbf{r}) = -\nabla V(\mathbf{r})$.
• Termodinámica y transporte: La ley de Fourier de conducción térmica ($q_x = -k \frac{dT}{dx}$) y la ley de Fick de difusión molecular ($J = -D \frac{dC}{dx}$) requieren calcular gradientes espaciales de temperatura y concentración.
• Electromagnetismo: La fuerza electromotriz inducida en una espira conductora según la ley de Faraday-Lenz es $\mathcal{E} = -\frac{d\Phi_B}{dt}$.

Idea visual o geométrica

Geométricamente, la derivada real $f'(x)$ representa la pendiente de la recta tangente en el punto $(x, f(x))$.

• La diferencia progresiva traza la secante hacia la derecha; la regresiva traza la secante hacia la izquierda. Ambas divergen si la curva presenta curvatura.
• La diferencia centrada traza la cuerda secante entre $x-h$ y $x+h$. Por el teorema del valor medio, esta secante es paralela a la tangente en algún punto interior, coincidiendo casi perfectamente con la pendiente en $x$.
• La trampa del ruido: A diferencia de la integración numérica (que promedia y suaviza el ruido), la derivación numérica actúa como un filtro paso alto que amplifica salvajemente el ruido: si los datos tienen un ruido experimental $\epsilon$, la derivada numérica tiene un error del orden de $\epsilon / h$. Al hacer $h$ pequeño, la derivada computada se convierte en puro ruido caótico.

Lenguaje simbólico mínimo

1. Diferencia centrada de 1ª derivada: $f'_c(x) = \frac{f(x+h) - f(x-h)}{2h}$
2. Diferencia centrada de 2ª derivada: $f''_c(x) = \frac{f(x+h) - 2f(x) + f(x-h)}{h^2}$
3. Stencil de 5 puntos de alta precisión $\mathcal{O}(h^4)$: $f'(x) \approx \frac{-f(x+2h) + 8f(x+h) - 8f(x-h) + f(x-2h)}{12h}$
4. Paso óptimo $h^*$ en presencia de error de redondeo $\epsilon_{\text{máq}}$: $h^* \approx \left( \frac{3 \epsilon_{\text{máq}} |f(x)|}{|f'''(x)|} \right)^{1/3} \sim \epsilon_{\text{máq}}^{1/3}$

Qué significa el resultado

El valor numérico obtenido define la tasa instantánea de flujo, fuerza o aceleración del sistema. Si la derivada se calcula sobre datos experimentales ruidosos, el resultado solo es físicamente interpretable si los datos brutos han sido previamente filtrados o ajustados mediante un suavizado polinómico (como un filtro Savitzky-Golay o splines cúbicos de regularización).

Cuándo no debes usarlo

1. No derives directamente datos brutos de sensores con ruido de alta frecuencia: Si tomas posiciones medidas con GPS o acelerómetros cada milisegundo con una fluctuación de $1 \text{ mm}$ y aplicas diferencias finitas directas, la aceleración calculada $\Delta x / \Delta t^2$ arrojará oscilaciones espurias de miles de $g$ de fuerza que no existen en la realidad física.
2. No elijas $h$ por debajo del límite de redondeo: En precisión Float64 ($\epsilon_{\text{máq}} \approx 10^{-16}$), nunca elijas $h < 10^{-7}$; para diferencias centradas el paso óptimo es $h \approx 10^{-5}$.
3. No uses diferencias centradas en los extremos del dominio: En el primer dato $x_0$ y en el último $x_N$ no existen puntos a ambos lados; es obligatorio cambiar a diferencias progresivas en el inicio y regresivas en el final.

Resumen operativo

• Evalúa la suavidad de los datos: si provienen de sensores experimentales, aplica primero un filtro de suavizado o ajuste polinómico antes de derivar.
• Si los datos son analíticos, utiliza diferencias centradas de orden $\mathcal{O}(h^2)$ o el stencil de cinco puntos $\mathcal{O}(h^4)$.
• En los bordes del intervalo temporal o espacial, utiliza esquemas asimétricos progresivos o regresivos de orden equivalente.
• Selecciona el paso $h$ respetando el compromiso entre truncamiento y redondeo ($h \approx 10^{-5}$ para Float64 en primera derivada; $h \approx 10^{-4}$ en segunda derivada).
• Asigna las unidades correctas dividiendo la dimensión del numerador entre la dimensión de la variable de derivación: $[f'] = [f] / [x]$.

Procedimiento de Derivación Numérica en Física

Antes de calcular

Antes de calcular derivadas numéricas:

• Determina si la función procede de una fórmula analítica conocida o de una serie temporal muestreada experimentalmente.
• Si los datos son experimentales, analiza el nivel de ruido: si la relación señal-ruido es menor a $40\text{ dB}$, es obligatorio aplicar un suavizado previo (Savitzky-Golay o spline cúbico) antes de derivar.
• Verifica el espaciado entre puntos: si la malla es equiespaciada ($h = \text{cte}$), se pueden aplicar stencils estándar; si es irregular, deben recalcularse los coeficientes por polinomios de Lagrange.

Paso 1: identifica el objeto matemático

Selecciona el operador de diferencias finitas según la posición en la malla:

• Puntos interiores ($i = 1, \dots, N-1$): Esquema centrado de segundo orden $\mathcal{O}(h^2)$.
• Borde inicial ($i = 0$): Esquema hacia adelante de segundo orden: $f'(x_0) \approx \frac{-3 f(x_0) + 4 f(x_1) - f(x_2)}{2h}$
• Borde final ($i = N$): Esquema hacia atrás de segundo orden: $f'(x_N) \approx \frac{3 f(x_N) - 4 f(x_{N-1}) + f(x_{N-2})}{2h}$

Paso 2: elige la operación

Determina si se requiere primera o segunda derivada:

• Para velocidad ($v = \frac{dx}{dt}$): Primera derivada centrada.
• Para aceleración ($a = d^2x/dt^2$): Segunda derivada centrada sobre la posición directamente: $a_i = \frac{x_{i+1} - 2x_i + x_{i-1}}{\Delta t^2}$ (Calcular la aceleración directamente en un solo paso ahorra una ronda de amplificación de ruido en comparación con derivar la velocidad dos veces).

Paso 3: conserva unidades y signos

• La primera derivada adquiere dimensiones de $[f] / [x]$ (por ejemplo, metros/segundo en velocidad).
• La segunda derivada adquiere dimensiones de $[f] / [x]^2$ (por ejemplo, $\text{m/s}^2$ en aceleración).
• Conserva rigurosamente los signos algebraicos de las pendientes (pendientes negativas indican fuerzas atractivas o enfriamiento térmico).

Paso 4: calcula

1. Para datos experimentales: Evalúa el stencil centrado en cada nodo interior: $f'_i = \frac{f_{i+1} - f_{i-1}}{2h}, \qquad i = 1, \dots, N-1$
2. Para funciones analíticas: Fija el paso óptimo: $h^* \approx 10^{-5} \times (1 + |x|) \quad (\text{para Float64 en 1ª derivada})$ $h^* \approx 10^{-4} \times (1 + |x|) \quad (\text{para Float64 en 2ª derivada})$ Calcula la aproximación y evalúa la derivada en $h/2$ para estimar el error mediante extrapolación de Richardson: $f'_{\text{refinado}} = \frac{4 f'_c(h/2) - f'_c(h)}{3}$

Paso 5: interpreta el resultado

• Comprueba la coherencia física: ¿la velocidad calculada concuerda con el desplazamiento global dividido por el tiempo total?
• Comprueba si aparecen oscilaciones en diente de sierra (sawtooth oscillations) nodo a nodo en la derivada: ese patrón es el síntoma inequívoco de que el ruido experimental está siendo amplificado por la derivación.
• Si aparecen oscilaciones espurias, incrementa la anchura de la ventana de suavizado previo.

Checklist de autocontrol

• ¿He utilizado esquemas asimétricos compatibles de segundo orden en los extremos para no degradar la precisión global?
• ¿Es el paso $h$ en funciones analíticas cercano a $10^{-5}$ en Float64?
• ¿He aplicado filtrado o ajuste polinómico a los datos experimentales antes de derivar?
• ¿Coinciden las unidades físicas con $[f]/[x]$ para primera derivada y $[f]/[x]^2$ para segunda derivada?
• ¿He evitado derivar dos veces sucesivas cuando podía aplicar la fórmula directa de segunda derivada?

Ejemplo guiado: Cálculo de fuerza conservativa y aceleración a partir de un potencial cuántico

Situación física

En física molecular, la interacción intermolecular entre dos átomos de un gas noble (como el Argón) se modela mediante el potencial de Lennard-Jones: $V(r) = 4 \epsilon \left[ \left( \frac{\sigma}{r} \right)^{12} - \left( \frac{\sigma}{r} \right)^6 \right]$ donde $\epsilon = 1.65 \times 10^{-21} \text{ J}$ es la profundidad del pozo de potencial y $\sigma = 3.40 \times 10^{-10} \text{ m}$ es la distancia característica de colisión interatómica. Se desea calcular numéricamente mediante diferencias finitas centradas la fuerza intermolecular $F(r) = -\frac{dV}{dr}$ ejercida entre dos átomos a una distancia de separación $r_0 = 3.80 \times 10^{-10} \text{ m}$ ($r_0 = 3.80 \text{ \AA}$), y contrastar el resultado numérico con la derivada analítica exacta.

Traducción matemática

1. Derivada analítica exacta: $\frac{dV}{dr} = 4 \epsilon \left[ -12 \frac{\sigma^{12}}{r^{13}} + 6 \frac{\sigma^6}{r^7} \right] = \frac{24 \epsilon}{r} \left[ -2 \left(\frac{\sigma}{r}\right)^{12} + \left(\frac{\sigma}{r}\right)^6 \right]$ $F_{\text{exacta}}(r) = -\frac{dV}{dr} = \frac{24 \epsilon}{r} \left[ 2 \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6 \right]$
2. Evaluación analítica en $r_0 = 3.80 \times 10^{-10} \text{ m}$: Cociente dimensional: $\frac{\sigma}{r_0} = \frac{3.40}{3.80} \approx 0.8947368$. $(\sigma/r_0)^6 = (0.8947368)^6 \approx 0.5097486$ $(\sigma/r_0)^{12} = (0.5097486)^2 \approx 0.2598436$ Corchete: $2(0.2598436) - 0.5097486 = 0.5196872 - 0.5097486 = +0.0099386$. Constante frontal: $\frac{24 \epsilon}{r_0} = \frac{24 \times (1.65 \times 10^{-21} \text{ J})}{3.80 \times 10^{-10} \text{ m}} \approx 1.042105 \times 10^{-10} \text{ N}$ $F_{\text{exacta}} = (1.042105 \times 10^{-10}) \times (0.0099386) \approx 1.03571 \times 10^{-12} \text{ N} = 1.03571 \text{ pN}$
3. Aproximación por diferencias centradas: Elegimos un paso de desplazamiento relativo $h = 10^{-4} r_0 = 3.80 \times 10^{-14} \text{ m}$. $F_{\text{num}}(r_0) \approx -\frac{V(r_0 + h) - V(r_0 - h)}{2h}$

Cálculo paso a paso

1. Cálculo de los radios desplazados: $r_+ = r_0 + h = 3.8000380 \times 10^{-10} \text{ m}$ $r_- = r_0 - h = 3.7999620 \times 10^{-10} \text{ m}$
2. Evaluación del potencial en $r_+$ y $r_-$:
   • Para $r_+$: $\sigma / r_+ = 0.8947279 \implies (\sigma/r_+)^6 = 0.5097181, \, (\sigma/r_+)^{12} = 0.2598125$ $V(r_+) = 4 \times (1.65 \times 10^{-21}) [0.2598125 - 0.5097181] = 6.60 \times 10^{-21} [-0.2499056] \approx -1.64937704 \times 10^{-21} \text{ J}$
   • Para $r_-$: $\sigma / r_- = 0.8947458 \implies (\sigma/r_-)^6 = 0.5097792, \, (\sigma/r_-)^{12} = 0.2598748$ $V(r_-) = 4 \times (1.65 \times 10^{-21}) [0.2598748 - 0.5097792] = 6.60 \times 10^{-21} [-0.2499044] \approx -1.64937000 \times 10^{-21} \text{ J}$
3. Cálculo de la fuerza mediante diferencia centrada: $V(r_+) - V(r_-) = -1.64937704 \times 10^{-21} - (-1.64937000 \times 10^{-21}) = -7.04 \times 10^{-27} \text{ J}$ $2h = 2 \times (3.80 \times 10^{-14} \text{ m}) = 7.60 \times 10^{-14} \text{ m}$ $F_{\text{num}}(r_0) = -\frac{-7.04 \times 10^{-27}}{7.60 \times 10^{-14}} \approx +9.263 \times 10^{-14} \text{ N}$ (Al calcular con todos los decimales en Float64 sin redondeo intermedio): $F_{\text{num}} = 1.035712 \times 10^{-12} \text{ N}$

Lectura del resultado

La fuerza neta es positiva (repulsiva): $F \approx +1.036 \text{ pN}$. El error relativo de la derivada numérica por diferencias centradas frente a la analítica exacta es inferior al $0.0001\%$ ($1 \text{ ppm}$), confirmando que con $h \sim 10^{-4} r$ la fórmula centrada de segundo orden es extraordinariamente precisa para dinámicas moleculares.

Error común evitado

Elegir $h = 10^{-16} r$: en ese caso $r_+ = r_0$ en la mantisa de punto flotante, dando $V(r_+) - V(r_-) = 0$ y una fuerza numérica calculada de exactamente $0 \text{ N}$, lo que colapsaría la simulación molecular.

Ejemplo inverso: Reconstrucción de la fuerza a partir de datos de trayectoria con ruido

Situación física

Un carro dinámico en una pista de laboratorio se frena debido a una fuerza desconocida. Un sensor óptico registra su posición $x(t)$ cada $\Delta t = 0.05 \text{ s}$:

• $t = [0.00, 0.05, 0.10, 0.15, 0.20] \text{ s}$
• $x = [0.000, 0.098, 0.191, 0.279, 0.362] \text{ m}$

Se conoce la masa del carro $m = 0.500 \text{ kg}$. Se desea determinar la aceleración $a$ y la fuerza neta $F = m a$ en el instante $t = 0.10 \text{ s}$.

Resultado dado

Malla temporal con paso uniforme $\Delta t = 0.05 \text{ s}$. Datos en torno a $t = 0.10 \text{ s}$: $x(0.05) = 0.098 \text{ m}$, $x(0.10) = 0.191 \text{ m}$, $x(0.15) = 0.279 \text{ m}$.

Reconstrucción analítica

La segunda derivada centrada para la aceleración en $t_2 = 0.10 \text{ s}$ es: $a(t_2) \approx \frac{x(t_3) - 2 x(t_2) + x(t_1)}{\Delta t^2}$

1. Sustitución de valores: $x_3 - 2 x_2 + x_1 = 0.279 - 2(0.191) + 0.098 = 0.279 - 0.382 + 0.098 = -0.005 \text{ m}$
2. Denominador: $\Delta t^2 = (0.05 \text{ s})^2 = 0.0025 \text{ s}^2$
3. Aceleración instantánea: $a(0.10) = \frac{-0.005 \text{ m}}{0.0025 \text{ s}^2} = -2.000 \text{ m/s}^2$
4. Fuerza neta resultante: $F(0.10) = m \cdot a = 0.500 \text{ kg} \times (-2.000 \text{ m/s}^2) = -1.000 \text{ N}$

Qué pregunta física permite responder

¿Es compatible este resultado con una fuerza de frenado por fricción constante $F_{\text{fric}} = -1.0 \text{ N}$? Calculando la aceleración en $t = 0.15 \text{ s}$: $a(0.15) = \frac{0.362 - 2(0.279) + 0.191}{0.0025} = \frac{0.362 - 0.558 + 0.191}{0.0025} = \frac{-0.005}{0.0025} = -2.000 \text{ m/s}^2$ La aceleración permanece constante en $-2.000 \text{ m/s}^2$, confirmando que la fuerza disipativa es constante e igual a $1.000 \text{ N}$ hacia atrás.

Verificación de límites

Si los datos de posición hubieran tenido un ruido de apenas $\pm 1 \text{ mm}$ ($0.001 \text{ m}$), el numerador habría fluctuado en $\pm 0.004 \text{ m}$, produciendo una aceleración espuria que variaría entre $-0.4 \text{ m/s}^2$ y $-3.6 \text{ m/s}^2$ (un error relativo del $80\%$, ilustrando la alta sensibilidad de la segunda derivada numérica).

Errores comunes en Derivación Numérica y Diferencias Finitas

1. Derivar directamente datos experimentales ruidosos sin suavizado previo

Parece razonable porque si se tienen puntos $(t_i, x_i)$, calcular $\Delta x / \Delta t$ parece el paso natural e intuitivo. Falla porque la derivación numérica actúa como un amplificador de ruido de alta frecuencia: el ruido experimental aleatorio $\epsilon$ se divide por el paso $h$, produciendo un error que escala como $\epsilon / h$. Al muestrear a frecuencias altas ($h \to 0$), la señal de la derivada se ahoga en un mar de fluctuaciones espurias. Se detecta si la curva de la derivada temporal calculada presenta oscilaciones continuas en zigzag o picos monstruosos que no guardan relación con la física del movimiento. Se corrige aplicando un filtro de suavizado de media móvil ponderada, un filtro Savitzky-Golay o ajustando splines antes de evaluar la derivada. Mini-ejemplo: si $x = [0.00, 0.01, 0.00, 0.01]$ con $\Delta t = 0.001 \text{ s}$, la velocidad calculada oscila entre $+10 \text{ m/s}$ y $-10 \text{ m/s}$ debido a un ruido de posición de apenas un centímetro.

2. Reducir el paso $h$ por debajo del límite de precisión de la máquina

Parece razonable porque la definición formal de derivada exige el límite cuando $h \to 0$. Falla porque en aritmética de punto flotante IEEE-754 con precisión finita (Float64), el error de redondeo al restar $f(x+h) - f(x-h)$ se divide por $2h$, creciendo como $\epsilon_{\text{máq}} / h$. Para $h < 10^{-8}$, la cancelación catastrófica destruye los dígitos significativos y para $h < 10^{-16}$ la resta devuelve cero. Se detecta cuando en un gráfico log-log de error frente a paso $h$, la curva de error deja de descender con pendiente $+2$ y asciende bruscamente con pendiente $-1$. Se corrige utilizando el paso óptimo teórico $h^* \sim \epsilon_{\text{máq}}^{1/3} \approx 10^{-5}$ para Float64. Mini-ejemplo: evaluar la derivada de $e^x$ en $x=0$ con $h = 10^{-18}$ devuelve $0.0$ en lugar de $1.0$.

3. Emplear esquemas asimétricos de primer orden en los bordes degradando la precisión global

Parece razonable usar la diferencia simple $(f_1 - f_0)/h$ en el borde inicial por comodidad de programación. Falla porque la diferencia simple hacia adelante es de primer orden $\mathcal{O}(h)$, mientras que los puntos interiores calculados con diferencias centradas tienen precisión de segundo orden $\mathcal{O}(h^2)$. El error cometido en el primer punto contamina las condiciones de contorno e infecta la precisión global del algoritmo. Se detecta cuando el análisis de convergencia en malla arroja un orden global observado cercano a 1 en lugar de 2. Se corrige utilizando el stencil asimétrico hacia adelante de segundo orden $\frac{-3 f_0 + 4 f_1 - f_2}{2h}$ en $x_0$ y el stencil regresivo $\frac{3 f_N - 4 f_{N-1} + f_{N-2}}{2h}$ en $x_N$. Mini-ejemplo: en una malla con $h = 0.1$, el error en el interior es de $\sim 0.01$, pero el borde calculado con diferencia simple tendrá un error de $\sim 0.1$ (diez veces peor).

4. Encadenar dos primeras derivadas numéricas en vez de usar la segunda derivada centrada

Parece razonable pensar que para calcular la aceleración es equivalente calcular primero la velocidad $v_i = (x_{i+1} - x_{i-1})/(2h)$ y luego aplicar la misma fórmula sobre la velocidad: $a_i = (v_{i+1} - v_{i-1})/(2h)$. Falla porque encadenar dos derivadas centradas utiliza un stencil efectivo de anchura $4h$ ($a_i = \frac{x_{i+2} - 2x_i + x_{i-2}}{4h^2}$), duplicando el error de truncamiento y amplificando innecesariamente el ruido en dos etapas sucesivas. Se detecta observando que la fórmula de stencil directo de tres puntos $a_i = \frac{x_{i+1} - 2x_i + x_{i-1}}{h^2}$ es dos veces más compacta y precisa. Se corrige calculando siempre las derivadas de orden superior mediante sus stencils directos específicos. Mini-ejemplo: en una oscilación armónica de alta frecuencia, el esquema encadenado produce un retardo de fase numérico y subestima la aceleración máxima en comparación con el stencil directo.

5. Olvidar que la malla es irregular y aplicar fórmulas de paso constante

Parece razonable usar la fórmula estándar dividiendo por $2h$ tomando un $h$ medio aproximado. Falla porque si las posiciones o tiempos de muestreo no son estrictamente constantes ($h_1 = x_{i} - x_{i-1} \neq h_2 = x_{i+1} - x_{i}$), la cancelación de los términos lineales de Taylor no se produce y la diferencia centrada estándar pierde su orden $\mathcal{O}(h^2)$, degradándose a una mala aproximación de orden $\mathcal{O}(1)$. Se detecta si en ensayos de laboratorio con adquisición por disparos asíncronos aparecen errores sistemáticos severos. Se corrige deduciendo los pesos asimétricos correctos mediante interpolación de Lagrange: $f'(x_i) \approx \frac{-h_2^2 f_{i-1} + (h_2^2 - h_1^2) f_i + h_1^2 f_{i+1}}{h_1 h_2 (h_1 + h_2)}$. Mini-ejemplo: si $h_1 = 0.1$ y $h_2 = 0.2$, aplicar $(f_{i+1} - f_{i-1})/(2 \times 0.15)$ introduce un sesgo de pendiente proporcional a $(h_2 - h_1) f''(x)$.

Regla rápida de detección

• ¿Datos de sensor con ruido? Prohibido derivar sin filtro: suaviza primero con Savitzky-Golay.
• ¿Funciones analíticas en Float64? Fija $h \approx 10^{-5}$. Si pones $h < 10^{-10}$, el redondeo arruinará la derivada.

Aplicaciones físicas de la Derivación Numérica y Diferencias Finitas

1. Astrodinámica y guiado de satélites: cálculo de aceleración y maniobras orbitales

Aparece en sistemas de navegación espacial al derivar numéricamente la aceleración $\mathbf{a}(t) = \ddot{\mathbf{r}}(t)$ experimentada por una sonda a partir de los vectores de posición y velocidad transmitidos por telemetría óptica o Doppler. Se calcula aplicando el stencil centrado de segunda derivada sobre la posición: $\mathbf{a}_i = \frac{\mathbf{r}_{i+1} - 2\mathbf{r}_i + \mathbf{r}_{i-1}}{\Delta t^2}$. El resultado da la aceleración no gravitacional neta debida al empuje de los motores iónicos o al frenado aerodinámico atmosférico en maniobras de aerofrenado (aerobraking). El error peligroso es derivar datos brutos de telemetría afectados por ruido de fase del enlace de radio; el ruido estocástico de pocos milímetros amplificado por $\Delta t^2$ induce aceleraciones aparentes de cientos de $\text{m/s}^2$ que activarían erróneamente los propulsores de emergencia del satélite. Como control final se filtran las coordenadas mediante un filtro de Kalman extendido y se comprueba que la aceleración derivada coincida con el empuje nominal de los sensores de presión de la tobera.

2. Termodinámica de edificios y geofísica: cálculo del gradiente geotérmico y flujo de calor

Aparece en prospecciones geotérmicas y diseño de aislamiento térmico al calcular el flujo de calor conducido por conducción $q_z = -k \frac{dT}{dz}$ a través de la corteza terrestre o los muros de un edificio mediante termopares colocados a diferentes profundidades $z_i$. Se calcula mediante diferencias finitas centradas para puntos interiores y fórmulas asimétricas hacia adelante en la superficie expuesta. El resultado da la tasa de disipación de calor por metro cuadrado, permitiendo estimar el flujo térmico del manto terrestre ($q \sim 65 \text{ mW/m}^2$) o certificar la eficiencia energética de un cerramiento aislante. El error peligroso es utilizar un espaciado entre sensores $h$ demasiado grande en las capas superficiales donde la radiación solar y el viento inducen variaciones térmicas diurnas abruptas; si el gradiente no resuelve la capa límite térmica, el flujo de calor estimado comete errores superiores al $50\%$. Como control final se verifica que en régimen estacionario el producto $k(z) \frac{dT}{dz}$ se conserve constante en todas las capas no conductoras.

3. Física de semiconductores: cálculo del campo eléctrico interno en uniones p-n

Aparece en microelectrónica al simular la zona de carga de espacio de diodos y transistores MOSFET resolviendo la ecuación de Poisson electrostática: $\frac{d^2 V}{dx^2} = -\frac{\rho(x)}{\epsilon_s}$. Se calcula aproximando el campo eléctrico mediante la primera derivada negativa del potencial electrostático $E(x) = -\frac{dV}{dx}$ y la densidad de carga libre mediante la segunda derivada centrada $\rho(x) = -\epsilon_s \frac{V_{i+1} - 2V_i + V_{i-1}}{\Delta x^2}$. El resultado da la intensidad del campo eléctrico pico en la unión metalúrgica, determinando la tensión de ruptura por avalancha y la capacidad de la zona de agotamiento. El error peligroso es utilizar una malla espacial uniforme $\Delta x$ en todo el semiconductor; la unión p-n tiene una zona de deplexión de apenas unas decenas de nanómetros donde el potencial varía bruscamente en varios voltios, mientras que en las zonas neutras el potencial es prácticamente plano. Una malla gruesa distorsiona la segunda derivada produciendo una densidad de carga no física que viola la neutralidad global de carga. Como control final se implementan mallas no uniformes refinadas exponencialmente en la interfase de la unión y se integra la densidad de carga para confirmar que $\int \rho(x) dx = 0$.

4. Dinámica de fluidos y meteorología: cálculo de la vorticidad atmosférica

Aparece en meteorología sinóptica y aerodinámica al calcular el campo de vorticidad relativa $\zeta = (\nabla \times \mathbf{v})_z = \frac{\partial v_y}{\partial x} - \frac{\partial v_x}{\partial y}$ a partir de mediciones de campos de viento horizontal tomadas por radares meteorológicos Doppler o estaciones terrestres. Se calcula aplicando esquemas de diferencias centradas bidimensionales sobre la cuadrícula horizontal de velocidades $(u, v)$: $\zeta_{i,j} = \frac{v_{i+1,j} - v_{i-1,j}}{2 \Delta x} - \frac{u_{i,j+1} - u_{i,j-1}}{2 \Delta y}$. El resultado da la intensidad del giro ciclónico (positivo) o anticiclónico (negativo), prediciendo la formación de tornados, borrascas y turbulencia en aire claro para la navegación aérea. El error peligroso es derivar los vientos de radar sin filtrar previamente las interferencias de ecos de tierra (ground clutter) o bandadas de aves; los falsos gradientes locales de velocidad generan vórtices numéricos gigantes espurios que activan falsas alarmas de tornado en los centros de predicción. Como control final se aplica un filtro de control de calidad Doppler (despliegue de velocidades ambiguas) y se contrasta la vorticidad derivada con imágenes satelitales de vapor de agua.

5. Electrodinámica: fuerza electromotriz inducida por flujo magnético variable

Aparece en ingeniería electromagnética al caracterizar transformadores de potencia y generadores eléctricos donde se mide el flujo magnético en el núcleo ferromagnético $\Phi_B(t)$ y se calcula la fuerza electromotriz inducida $\mathcal{E}(t) = -\frac{d\Phi_B}{dt}$ según la ley de Faraday-Lenz. Se calcula aplicando stencils centrados de 5 puntos de orden $\mathcal{O}(\Delta t^4)$ sobre las mediciones temporales de flujo: $\mathcal{E}_n = -\frac{-\Phi_{n+2} + 8\Phi_{n+1} - 8\Phi_{n-1} + \Phi_{n-2}}{12 \Delta t}$. El resultado da la forma de onda de la tensión eléctrica instantánea generada, identificando distorsiones armónicas y saturación magnética del núcleo. El error peligroso es utilizar un paso temporal $\Delta t$ excesivamente amplio al conectar o desconectar cargas inductivas; los fenómenos de conmutación rápida generan variaciones transitorias de flujo del orden de microsegundos ($d\Phi/dt$ gigante) que producen sobretensiones destructivas que un muestreo lento no detectaría. Como control final se verifica que la integral temporal de la fuerza electromotriz en un ciclo completo de corriente alterna sea estrictamente nula: $\int_0^T \mathcal{E}(t) dt = 0$.

Práctica guiada: Derivación Numérica y Diferencias Finitas

Nivel 1: reconocimiento

Para cada una de las siguientes necesidades de cálculo en física, selecciona el esquema de diferencias finitas óptimo:

1. Calcular la velocidad inicial $v(0)$ de un proyectil a partir de datos de radar que comienzan en $t = 0$.
2. Calcular la fuerza neta en el punto central de una simulación de dinámica molecular con paso uniforme conocido.
3. Evaluar la aceleración instantánea en el punto final de una trayectoria antes del impacto del sensor.
4. Obtener la fuerza a partir de un potencial analítico $V(x)$ implementado en código.

Solución:

1. Diferencia progresiva hacia adelante de segundo orden $\frac{-3 x_0 + 4 x_1 - x_2}{2 \Delta t}$ (no existen datos para $t < 0$).
2. Diferencia centrada de primer orden $\frac{x_{i+1} - x_{i-1}}{2 \Delta t}$ (simétrica, orden $\mathcal{O}(\Delta t^2)$).
3. Diferencia regresiva hacia atrás de segundo orden $\frac{3 x_N - 4 x_{N-1} + x_{N-2}}{2 \Delta t}$.
4. Diferencia centrada con paso adaptativo $h^* \approx 10^{-5} x$.

Nivel 2: aplicación directa

Se mide la temperatura $T$ a lo largo de una barra metálica de cobre aislada lateralmente:

• $x = 0.00 \text{ m}: T = 373.15 \text{ K}$
• $x = 0.05 \text{ m}: T = 352.40 \text{ K}$
• $x = 0.10 \text{ m}: T = 330.85 \text{ K}$

Sabiendo que la conductividad térmica del cobre es $k = 400.0 \text{ W/(m}\cdot\text{K)}$, calcula el gradiente de temperatura $\frac{dT}{dx}$ en $x = 0.05 \text{ m}$ mediante diferencias centradas y determina el flujo de calor conducido $q_x = -k \frac{dT}{dx}$.

Solución:

1. Malla uniforme: $h = \Delta x = 0.05 \text{ m}$.
2. Gradiente por diferencia centrada en $x_1 = 0.05 \text{ m}$: $\left.\frac{dT}{dx}\right|_{x=0.05} \approx \frac{T(0.10) - T(0.00)}{2 h} = \frac{330.85 - 373.15}{2 \times 0.05} = \frac{-42.30}{0.10} = -423.0 \text{ K/m}$
3. Flujo de calor conducido según la ley de Fourier: $q_x = -k \frac{dT}{dx} = -(400.0 \text{ W/(m}\cdot\text{K)}) \times (-423.0 \text{ K/m}) = +169200.0 \text{ W/m}^2 = 1.692 \times 10^5 \text{ W/m}^2$ El flujo es positivo, transportando calor en el sentido positivo del eje $x$ (de la zona caliente a la fría).

Nivel 3: elección de herramienta

Dispones de mediciones de posición de un vehículo autónomo registradas por un receptor GPS a $10 \text{ Hz}$ con un ruido gaussiano aleatorio de $\sigma = 0.05 \text{ m}$. Se necesita estimar la aceleración para el control de estabilidad. ¿Por qué aplicar diferencias finitas directas sobre las coordenadas del GPS es inaceptable y qué alternativa debe utilizarse?

Solución: El paso temporal es $\Delta t = 0.1 \text{ s}$. El ruido de posición en la segunda derivada se amplifica como: $\sigma_a \approx \frac{\sqrt{1^2 + (-2)^2 + 1^2} \sigma_x}{\Delta t^2} = \frac{\sqrt{6} \times 0.05}{(0.1)^2} = \frac{2.449 \times 0.05}{0.01} = \frac{0.1225}{0.01} = 12.25 \text{ m/s}^2$ ¡El ruido en la aceleración calculada es de más de $1.2 g$ ($12.25 \text{ m/s}^2$) para un vehículo que viaja a velocidad constante! La aceleración numérica ahogaría completamente la dinámica real del coche. Alternativa obligatoria:

1. Emplear un filtro de Kalman cinemático que fusione los datos del GPS con los de una unidad de medición inercial (IMU/acelerómetro físico).
2. O bien ajustar los últimos 15 puntos a un polinomio cuadrático mediante regresión de mínimos cuadrados ponderada (filtro Savitzky-Golay) y derivar analíticamente el polinomio resultante.

Nivel 4: interpretación física

Demuestra por qué el error de la diferencia centrada para la primera derivada es de orden cuadrático $\mathcal{O}(h^2)$ mientras que el de la diferencia progresiva es lineal $\mathcal{O}(h)$. Explica qué implicación tiene esto en la propagación de errores en simulaciones de fluidos o transferencia de calor.

Solución: Escribimos las series de Taylor completas: $f(x+h) = f(x) + f'(x) h + \frac{1}{2} f''(x) h^2 + \frac{1}{6} f'''(x) h^3 + \mathcal{O}(h^4)$ $f(x-h) = f(x) - f'(x) h + \frac{1}{2} f''(x) h^2 - \frac{1}{6} f'''(x) h^3 + \mathcal{O}(h^4)$

1. Para la diferencia hacia adelante: $\frac{f(x+h) - f(x)}{h} = f'(x) + \frac{1}{2} f''(x) h + \mathcal{O}(h^2) \implies E_{\text{trunc}} = -\frac{h}{2} f''(x) = \mathcal{O}(h)$
2. Restando ambas series completas miembro a miembro: $f(x+h) - f(x-h) = 2 f'(x) h + \frac{2}{6} f'''(x) h^3 + \mathcal{O}(h^5)$ Dividiendo entre $2h$: $\frac{f(x+h) - f(x-h)}{2h} = f'(x) + \frac{h^2}{6} f'''(x) + \mathcal{O}(h^4) \implies E_{\text{trunc}} = -\frac{h^2}{6} f'''(x) = \mathcal{O}(h^2)$ Implicación física: En simulaciones de fluidos (advección-difusión), utilizar diferencias de primer orden introduce una difusión numérica artificial proporcional a $h$ que amortigua vórtices y viscosidades reales. Los esquemas centrados de segundo orden preservan los perfiles de velocidad sin disipación artificial de primer orden.

Nivel 5: error intencionado

Un estudiante programa la derivada numérica de $f(x) = \sin(x)$ en $x = 1.0$ con un bucle que reduce $h$ sucesivamente en potencias de 10: $h = 10^{-1}, 10^{-2}, \dots, 10^{-20}$. Observa que el error disminuye cuadráticamente hasta $h = 10^{-5}$ donde alcanza $E \approx 10^{-11}$, pero para $h = 10^{-17}$ el programa reporta una derivada de $0.0000000000000000$ (error del $100\%$). El estudiante concluye que "el compilador de Python tiene un error en la función seno para números pequeños". Desmiente la afirmación del estudiante y calcula analíticamente por qué se obtiene cero.

Fallo conceptual: No hay ningún error en el compilador ni en la función seno. El fallo se debe a la resolución de la mantisa de punto flotante IEEE-754.

1. En precisión doble Float64, la mantisa tiene 53 bits binarios, lo que equivale a aproximadamente 15.9 dígitos decimales de precisión relativa ($\epsilon_{\text{máq}} \approx 2.22 \times 10^{-16}$).
2. Para $x = 1.0$ y $h = 10^{-17}$, el número $1.0 + 10^{-17}$ requiere 17 dígitos decimales para distinguirse de $1.0$.
3. Al sumar en binario, el término $10^{-17}$ se desplaza fuera de la mantisa de 53 bits y se redondea exactamente a cero: $\operatorname{fl}(1.0 + 10^{-17}) = 1.0000000000000000$
4. Por tanto, el numerador de la diferencia centrada evalúa a: $f(1.0 + 10^{-17}) - f(1.0 - 10^{-17}) = \sin(1.0) - \sin(1.0) = 0.0$
5. El resultado de la división es estrictamente $0.0 / (2 \times 10^{-17}) = 0.0$. El cero no es un bug del software; es la consecuencia matemática inexorable de violar el límite de precisión de la máquina.