Integración numérica

Cuadratura numérica mediante reglas del trapecio, Simpson y métodos de orden superior para integrales.

Qué es

La integración numérica, históricamente denominada cuadratura numérica, es el conjunto de métodos algorítmicos que permite calcular el valor numérico aproximado de una integral definida: $I = \int_a^b f(x) dx$ cuando la función integrando $f(x)$ carece de primitiva analítica elemental (como $e^{-x^2}$, $\frac{\sin x}{x}$, o $\sqrt{1 - k^2 \sin^2\theta}$), o cuando $f(x)$ se conoce únicamente como un conjunto discreto de puntos muestreados experimentalmente en el laboratorio.

Los esquemas de cuadratura se estructuran en dos grandes categorías:

1. Fórmulas compuestas de Newton-Cotes (mallas equiespaciadas con paso $h = \frac{b-a}{N}$):
   • Regla del Trapecio compuesta: Conecta los puntos mediante segmentos rectilíneos. Su error de truncamiento global es de segundo orden: $I_{\text{trap}} = \frac{h}{2} \left[ f(a) + 2 \sum_{i=1}^{N-1} f(x_i) + f(b) \right], \qquad E_{\text{trunc}} = -\frac{(b-a) h^2}{12} f''(\xi) = \mathcal{O}(h^2)$
   • Regla de Simpson 1/3 compuesta ($N$ par): Interpola arcos parabólicos cuadráticos sobre pares de subintervalos contiguos. Su error global es de cuarto orden: $I_{\text{Simp}} = \frac{h}{3} \left[ f(x_0) + 4 \sum_{i=1,3,\dots}^{N-1} f(x_i) + 2 \sum_{i=2,4,\dots}^{N-2} f(x_i) + f(x_N) \right], \qquad E_{\text{trunc}} = -\frac{(b-a) h^4}{180} f^{(4)}(\xi) = \mathcal{O}(h^4)$
   • Aceleración de Romberg: Aplica extrapolación de Richardson sobre la regla del trapecio con pasos $h, h/2, h/4, \dots$ eliminando sistemáticamente los términos de error $h^2, h^4, h^6$.
2. Cuadratura Gaussiana (nodos óptimos no equiespaciados):
   • Gauss-Legendre: Elige estratégicamente tanto las posiciones de los nodos $x_i$ como los pesos $w_i$ correspondientes a las raíces de los polinomios ortogonales de Legendre $P_n(x)$. Con apenas $n$ evaluaciones de la función, integra de forma exacta cualquier polinomio de grado hasta $2n - 1$: $\int_{-1}^1 f(t) dt \approx \sum_{i=1}^n w_i f(t_i)$

Por qué aparece en física

La integración está indisolublemente ligada a los principios acumulativos de conservación y campos en física:

• Mecánica clásica: El trabajo mecánico realizado por una fuerza variable no lineal es $W = \int_{x_1}^{x_2} F(x) dx$; el centro de masa de cuerpos continuos heterogéneos y los momentos de inercia requieren integrar distribuciones volumétricas de densidad $\rho(\mathbf{r})$.
• Electromagnetismo: El potencial electrostático y el campo magnético en el espacio circundante a distribuciones continuas de corriente se calculan mediante integrales de Coulomb y Biot-Savart, las cuales casi nunca tienen soluciones analíticas para geometrías reales de bobinas o electrodos.
• Termodinámica y física estadística: La ley de radiación de Stefan-Boltzmann ($u = a T^4$) surge integrando la densidad espectral de Planck en todas las frecuencias $\int_0^\infty \frac{x^3}{e^x - 1} dx = \frac{\pi^4}{15}$; las funciones de partición cuánticas y capacidades caloríficas de Debye involucran integrales numéricas continuas.
• Oscilaciones no lineales: El periodo exacto de oscilación de un péndulo simple con amplitudes angulares grandes $\theta_0$ es una integral elíptica completa de primera especie $T = 4 \sqrt{\frac{L}{g}} \int_0^{\pi/2} \frac{d\phi}{\sqrt{1 - \sin^2(\theta_0/2) \sin^2\phi}}$, imposible de resolver con funciones trigonométricas estándar.

Idea visual o geométrica

Geométricamente, la integral representa el área neta con signo bajo la curva $f(x)$ entre las cotas $a$ y $b$.

• En el Trapecio, el área real se cubre con un mosaico de trapecios contiguos cuyas tapas superiores forman una poligonal secante.
• En Simpson 1/3, la superficie se aproxima mediante una secuencia de franjas paraboloides que abrazan la curvatura local de la función mucho más de cerca.
• En la Cuadratura de Gauss, los puntos de evaluación se concentran hacia los extremos del intervalo donde la curvatura polinómica suele inducir mayor divergencia, maximizando la precisión global sin necesidad de densificar la malla global uniformemente.

Lenguaje simbólico mínimo

1. Regla del Trapecio compuesta: $T(h) = h \left[ \frac{f(a) + f(b)}{2} + \sum_{i=1}^{N-1} f(a + i h) \right]$
2. Regla de Simpson 1/3 compuesta: $S(h) = \frac{h}{3} \left[ f(x_0) + 4 \sum_{i=1,3,\dots}^{N-1} f(x_i) + 2 \sum_{i=2,4,\dots}^{N-2} f(x_i) + f(x_N) \right]$
3. Fórmula de Romberg de extrapolación: $R_{k,j} = \frac{4^j R_{k, j-1} - R_{k-1, j-1}}{4^j - 1}$
4. Transformación de coordenadas a $[-1, 1]$ para Gauss-Legendre: $x = \frac{b - a}{2} t + \frac{a + b}{2}, \qquad dx = \frac{b - a}{2} dt$ $\int_a^b f(x) dx = \frac{b - a}{2} \int_{-1}^1 f\left( \frac{b-a}{2}t + \frac{a+b}{2} \right) dt$

Qué significa el resultado

El valor escalar resultante expresa la magnitud física acumulada (julios de energía, voltios de potencial, segundos de periodo, culombios de carga). El orden de convergencia asintótica asegura que al duplicar el número de subintervalos $N$:

• En el Trapecio, el error se divide aproximadamente por $2^2 = 4$.
• En Simpson 1/3, el error se divide aproximadamente por $2^4 = 16$.
• Con Gauss-Legendre de orden $n=5$, se obtienen con frecuencia 12 o 14 dígitos significativos correctos con apenas 5 evaluaciones del integrando.

Cuándo no debes usarlo

1. Integrandos con singularidades integrables en los límites: Si $f(x) \to \infty$ en uno de los extremos (por ejemplo, $\int_0^1 \frac{1}{\sqrt{x}} dx$), evaluar Newton-Cotes cerrado en $x=0$ produce división por cero inmediata. Debe realizarse un cambio de variable analítico que regularice la singularidad o emplear cuadraturas abiertas de Gauss-Jacobi.
2. Integrandos altamente oscilatorios: Si $f(x) = g(x) \sin(\omega x)$ con frecuencia angular ultra-alta $\omega \gg 1$, las reglas estándar requieren mallas con millones de puntos para resolver cada ciclo de onda ($h \ll 2\pi/\omega$). En tales casos deben emplearse métodos de Filon o cuadratura asintótica.
3. Integrales multidimensionales de alta dimensión ($D > 4$): En mecánica estadística con decenas de partículas o integrales de camino de Feynman, la cuadratura en malla sufre la maldición de la dimensionalidad ($N^D$ evaluaciones). En esas dimensiones la integración por Monte Carlo o Monte Carlo cuántico es el único método computacionalmente viable.

Resumen operativo

• Examina el integrando $f(x)$: ¿es suave y diferenciable? ¿presenta discontinuidades o asíntotas en el intervalo $[a, b]$?
• Si dispones de datos experimentales tabulados a intervalos regulares, utiliza la regla de Simpson compuesta si $N$ es par, o la regla del trapecio si $N$ es impar o la malla es heterogénea.
• Si dispones de la expresión analítica continua y la función es suave, implementa Gauss-Legendre o cuadratura adaptativa de Simpson para máxima eficiencia.
• Aplica el test de duplicación de malla ($N$ y $2N$) para estimar el error residual a posteriori mediante extrapolación de Richardson.
• Multiplica el resultado numérico final por las unidades físicas compuestas correspondientes al producto $[f] \cdot [x]$.

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

Antes de calcular

Antes de calcular una integral numérica:

• Comprueba si el intervalo de integración $[a, b]$ contiene asíntotas verticales ($f(x) \to \pm \infty$) o discontinuidades de salto.
• Si los datos proceden de un experimento, verifica si el espaciado temporal o espacial $\Delta x$ es constante o variable.
• Si la función es analítica y tiene simetrías (par o impar), simplifica el dominio de integración antes de computar.

Paso 1: identifica el objeto matemático

Determina la técnica adecuada según la naturaleza del integrando:

• Datos tabulados equiespaciados con $N$ par: Regla de Simpson 1/3 compuesta.
• Datos tabulados con espaciado irregular: Regla del Trapecio no uniforme.
• Función matemática continua suave sin singularidades: Cuadratura de Gauss-Legendre de $n$ puntos.

Paso 2: elige la operación

Para la regla de Simpson 1/3 compuesta con $N$ subintervalos ($N$ par, $N+1$ puntos):

1. Calcula el paso de malla: $h = \frac{b - a}{N}$.
2. Genera los puntos de muestreo: $x_i = a + i h$ para $i = 0, 1, \dots, N$.
3. Agrupa los sumandos ponderados por coeficientes de Newton-Cotes ($1, 4, 2, 4, 2, \dots, 4, 1$).

Paso 3: conserva unidades y signos

• Los valores del integrando $f(x_i)$ tienen las dimensiones físicas de la función $[f]$.
• El paso de discretización $h$ tiene las dimensiones físicas de la variable de integración $[x]$.
• El resultado numérico compuesto tiene estrictamente unidades de $[f] \times [x]$ (por ejemplo, si $F$ está en newtons y $x$ en metros, el resultado tiene unidades exactas de $\text{N}\cdot\text{m} = \text{J}$).
• Respeta los signos algebraicos si los límites están invertidos: $\int_b^a f(x) dx = -\int_a^b f(x) dx$.

Paso 4: calcula

1. Evaluación de Simpson 1/3: $I(h) = \frac{h}{3} \left[ f(a) + 4 \sum_{k=1}^{N/2} f(x_{2k-1}) + 2 \sum_{k=1}^{N/2 - 1} f(x_{2k}) + f(b) \right]$
2. Evaluación con malla duplicada ($2N$ subintervalos): Calcula $I(h/2)$ utilizando el doble de puntos.
3. Estimación del error a posteriori (extrapolación de Richardson): $E_{\text{aprox}} \approx \frac{I(h/2) - I(h)}{15}$ El valor refinado de orden superior es: $I_{\text{refinado}} = I(h/2) + \frac{I(h/2) - I(h)}{15}$

Paso 5: interpreta el resultado

• Compara la cota de error $E_{\text{aprox}}$ con la tolerancia exigida por el problema físico (típicamente $10^{-6}$ en ingeniería o $10^{-12}$ en metrología).
• Si la estimación del error es mayor que la tolerancia, duplica nuevamente el número de puntos ($4N$) hasta alcanzar convergencia estable.
• Expresa el resultado final con sus correspondientes unidades físicas del SI.

Checklist de autocontrol

• ¿Es el número de subintervalos $N$ estrictamente par para aplicar la regla de Simpson 1/3?
• ¿He multiplicado el factor común por $h/3$ y no por $h/2$?
• ¿Coinciden las unidades del resultado con el producto $[f] \cdot [x]$?
• ¿He comprobado que el integrando no tenga una división por cero oculta dentro del intervalo $[a, b]$?
• ¿Disminuye el error en un factor cercano a 16 al duplicar los puntos en Simpson?

Ejemplo guiado: Periodo de oscilación de un péndulo no lineal para grandes amplitudes

Situación física

Un péndulo simple rígido de longitud $L = 1.000 \text{ m}$ en un campo gravitatorio de $g = 9.80665 \text{ m/s}^2$ se libera desde el reposo con un ángulo inicial de gran amplitud $\theta_0 = 60.0^\circ = \frac{\pi}{3} \text{ rad}$. La aproximación estándar de pequeñas oscilaciones ($\,\sin\theta \approx \theta$) predice: $T_0 = 2\pi \sqrt{\frac{L}{g}} = 2\pi \sqrt{\frac{1.000}{9.80665}} \approx 2.006409 \text{ s}$ Sin embargo, para $\theta_0 = 60^\circ$ la no linealidad del seno retrasa sustancialmente la oscilación. La conservación exacta de la energía mecánica conduce a la integral definida elíptica: $T = 4 \sqrt{\frac{L}{g}} \int_0^{\pi/2} \frac{d\phi}{\sqrt{1 - k^2 \sin^2\phi}}, \qquad k = \sin\left(\frac{\theta_0}{2}\right) = \sin(30^\circ) = 0.500$ Se desea calcular el periodo exacto $T$ mediante la regla de Simpson 1/3 compuesta con $N = 4$ subintervalos.

Traducción matemática

La constante frontal es: $C = 4 \sqrt{\frac{L}{g}} = 4 \sqrt{\frac{1.000}{9.80665}} = 4 \times 0.3193297 \approx 1.277319 \text{ s}$ El integrando a evaluar en $[0, \pi/2]$ es: $f(\phi) = \frac{1}{\sqrt{1 - (0.5)^2 \sin^2\phi}} = \frac{1}{\sqrt{1 - 0.25 \sin^2\phi}}$ Para $N = 4$ subintervalos, el paso de malla es: $h = \frac{\pi/2 - 0}{4} = \frac{\pi}{8} \approx 0.3926991 \text{ rad}$ Los 5 nodos de evaluación son $\phi = [0, \pi/8, 2\pi/8, 3\pi/8, 4\pi/8] = [0, \pi/8, \pi/4, 3\pi/8, \pi/2]$.

Cálculo paso a paso

1. Evaluación de la función en los nodos:
   • $\phi_0 = 0$: $\sin(0) = 0 \implies f(\phi_0) = \frac{1}{\sqrt{1 - 0}} = 1.0000000$
   • $\phi_1 = \pi/8 \approx 0.3926991$: $\sin(\pi/8) \approx 0.3826834 \implies \sin^2 = 0.1464466$ $f(\phi_1) = \frac{1}{\sqrt{1 - 0.25(0.1464466)}} = \frac{1}{\sqrt{1 - 0.0366117}} = \frac{1}{\sqrt{0.9633883}} \approx 1.0188048$
   • $\phi_2 = \pi/4 \approx 0.7853982$: $\sin(\pi/4) = \frac{\sqrt{2}}{2} \implies \sin^2 = 0.5000000$ $f(\phi_2) = \frac{1}{\sqrt{1 - 0.25(0.5000000)}} = \frac{1}{\sqrt{0.875}} \approx 1.0690450$
   • $\phi_3 = 3\pi/8 \approx 1.1780972$: $\sin(3\pi/8) \approx 0.9238795 \implies \sin^2 = 0.8535534$ $f(\phi_3) = \frac{1}{\sqrt{1 - 0.25(0.8535534)}} = \frac{1}{\sqrt{1 - 0.2133883}} = \frac{1}{\sqrt{0.7866117}} \approx 1.1274944$
   • $\phi_4 = \pi/2 \approx 1.5707963$: $\sin(\pi/2) = 1.0 \implies \sin^2 = 1.0000000$ $f(\phi_4) = \frac{1}{\sqrt{1 - 0.25(1.0)}} = \frac{1}{\sqrt{0.75}} \approx 1.1547005$
2. Aplicación de la regla de Simpson 1/3: $I_{\text{Simp}} = \frac{h}{3} \left[ f(\phi_0) + 4 f(\phi_1) + 2 f(\phi_2) + 4 f(\phi_3) + f(\phi_4) \right]$ Suma ponderada interna: $1.0000000 + 4(1.0188048) + 2(1.0690450) + 4(1.1274944) + 1.1547005$ $= 1.0000000 + 4.0752192 + 2.1380900 + 4.5099776 + 1.1547005 = 12.8779873$ Multiplicación por $h/3$: $I_{\text{Simp}} = \frac{0.3926991}{3} \times 12.8779873 \approx 0.1308997 \times 12.8779873 \approx 1.685725$
3. Cálculo del periodo final: $T = C \times I_{\text{Simp}} = 1.277319 \times 1.685725 \approx 2.15321 \text{ s}$

Lectura del resultado

El periodo exacto a $60^\circ$ es $T \approx 2.153 \text{ s}$, lo que representa un incremento del $+7.3\%$ respecto a la fórmula armónica lineal ($2.006 \text{ s}$). Con apenas $N=4$ intervalos, la regla de Simpson reproduce el valor exacto de la integral elíptica ($K(0.5) \approx 1.685750$) con un error relativo de apenas el $0.0015\%$.

Error común evitado

Utilizar la regla del trapecio con $N=4$: el trapecio daría un valor significativamente subestimado con un error casi 20 veces mayor que Simpson para el mismo número de puntos calculados.

Ejemplo inverso: Trabajo mecánico y reconstrucción del potencial a partir de mediciones de fuerza

Situación física

En un acelerador de iones se mide la fuerza electrostática no lineal de confinamiento $F(x)$ ejercida sobre una partícula cargada a lo largo del eje $x$ entre $x = 0.0 \text{ m}$ y $x = 0.6 \text{ m}$. Los sensores de fuerza proporcionan lecturas discretas espaciadas uniformemente cada $\Delta x = 0.1 \text{ m}$:

• $x = [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6] \text{ m}$
• $F = [0.00, 1.95, 3.68, 4.97, 5.64, 5.51, 4.32] \text{ N}$

Se conoce que el trabajo total realizado por la fuerza es $W = \int_0^{0.6} F(x) dx$.

Resultado dado

Malla discreta de 7 puntos ($N=6$ subintervalos pares, paso $h = 0.1 \text{ m}$).

Reconstrucción analítica

Como $N=6$ es estrictamente par, la regla de Simpson 1/3 compuesta es óptima: $W = \frac{h}{3} [F_0 + 4(F_1 + F_3 + F_5) + 2(F_2 + F_4) + F_6]$

1. Extremos: $F_0 + F_6 = 0.00 + 4.32 = 4.32 \text{ N}$
2. Nodos impares (peso 4): $F_1 + F_3 + F_5 = 1.95 + 4.97 + 5.51 = 12.43 \text{ N} \implies 4 \times 12.43 = 49.72 \text{ N}$
3. Nodos pares interiores (peso 2): $F_2 + F_4 = 3.68 + 5.64 = 9.32 \text{ N} \implies 2 \times 9.32 = 18.64 \text{ N}$
4. Suma total: $\sum = 4.32 + 49.72 + 18.64 = 72.68 \text{ N}$
5. Multiplicación: $W = \frac{0.1 \text{ m}}{3} \times 72.68 \text{ N} = \frac{7.268}{3} \approx 2.4227 \text{ J}$

Qué pregunta física permite responder

Si la partícula parte del reposo en $x=0$ con masa $m = 2.0 \times 10^{-3} \text{ kg}$, ¿cuál es su velocidad al alcanzar $x = 0.6 \text{ m}$? Por el teorema de las fuerzas vivas: $W = \Delta E_c = \frac{1}{2} m v_f^2 - 0 \implies v_f = \sqrt{\frac{2 W}{m}} = \sqrt{\frac{2 \times 2.4227}{2.0 \times 10^{-3}}} = \sqrt{2422.7} \approx 49.22 \text{ m/s}$

Verificación de límites

Si se hubiera empleado la regla trapezoidal ordinaria, el resultado habría sido $W_{\text{trap}} = 0.1 \times [0 + 1.95 + 3.68 + 4.97 + 5.64 + 5.51 + 2.16] = 2.391 \text{ J}$, subestimando la energía transferida en más del $1.3\%$ debido a la curvatura cóncava de la fuerza.

Errores comunes en Integración Numérica (Cuadratura)

1. Aplicar la regla de Simpson 1/3 a un número impar de subintervalos

Parece razonable porque uno asume que el software o fórmula simplemente toma los puntos y promedia. Falla porque la regla de Simpson 1/3 exige emparejar los subintervalos de dos en dos para ajustar una parábola en cada par ($x_0, x_1, x_2$), lo que exige estrictamente un número par de subintervalos $N$ (es decir, un número impar de puntos $N+1$). Si $N$ es impar, el último intervalo queda desemparejado. Se detecta si el algoritmo lanza un error de índice de matriz o si evalúa el último término con pesos incorrectos. Se corrige comprobando la paridad de $N$; si $N$ es impar, se aplica Simpson 1/3 en los primeros $N-1$ intervalos y la regla del trapecio (o Simpson 3/8) en el último intervalo. Mini-ejemplo: intentar aplicar Simpson 1/3 a 5 subintervalos ($N=5$, 6 puntos) deja un subintervalo sin parábola; forzar los coeficientes produce errores de truncamiento del $10\%$.

2. Evaluar reglas de Newton-Cotes cerradas en singularidades de frontera

Parece razonable porque muchas integrales físicas impropias convergen a valores finitos (como $\int_0^1 \frac{dx}{\sqrt{x}} = 2$). Falla porque las reglas cerradas (Trapecio y Simpson estándar) evalúan obligatoriamente el integrando en los extremos exactos $a$ y $b$. En $x=0$, $1/\sqrt{0}$ produce división por cero y detiene la ejecución. Se detecta mediante excepciones en tiempo de ejecución por división entre cero o valores infinitos. Se corrige mediante regularización analítica por sustitución ($x = u^2$), utilizando cuadraturas abiertas (como punto medio compuesto) o cuadraturas de Gauss que jamás evalúan los extremos de frontera. Mini-ejemplo: para $\int_0^1 x^{-1/2} dx$, el cambio $x = t^2$ produce $\int_0^1 2 dt = 2$, eliminando la singularidad por completo.

3. Discretizar integrales altamente oscilatorias con mallas demasiado gruesas

Parece razonable porque si el integrando oscila entre $-1$ y $+1$, uno podría pensar que las oscilaciones se cancelan aproximadamente. Falla porque si el paso de malla $h$ es mayor que el periodo de oscilación ($h > \lambda_{\text{onda}}$), se produce un falso muestreo (aliasing espacial) donde los nodos coinciden sistemáticamente con picos positivos o negativos, calculando un área gigante no física en vez de la cancelación real. Se detecta cuando al aumentar ligeramente el número de puntos la integral calculada cambia violentamente de signo o de magnitud. Se corrige exigiendo al menos 10 a 20 puntos por cada ciclo de oscilación completa ($h \leq \frac{2\pi}{10 \omega}$) o utilizando métodos analítico-numéricos de Filon. Mini-ejemplo: integrar $\int_0^{100\pi} \sin(10 x) dx$ con paso $h = \pi/5$ evalúa siempre $\sin(2\pi k) = 0$, arrojando $I=0$ en todo momento por pura casualidad aritmética, pero con $h = \pi/10 + 0.01$ arrojará valores caóticos descontrolados.

4. Confundir el orden de convergencia local de un subintervalo con el orden global

Parece razonable porque en la deducción elemental se muestra que el error en un único panel de Simpson es $\mathcal{O}(h^5)$. Falla porque para cubrir un intervalo fijo $[a, b]$ se deben sumar $N = (b-a)/h$ paneles. Al sumar $N$ términos de orden $\mathcal{O}(h^5)$, el error acumulado global se degrada en una potencia de $h$: $\sum_{i=1}^N \mathcal{O}(h^5) \approx N \mathcal{O}(h^5) = \frac{b-a}{h} \mathcal{O}(h^5) = \mathcal{O}(h^4)$. Se detecta cuando un analista asume que duplicar los puntos reducirá el error en $2^5 = 32$ veces en vez del factor real $2^4 = 16$. Se corrige utilizando siempre el orden global $\mathcal{O}(h^2)$ para el trapecio compuesto y $\mathcal{O}(h^4)$ para Simpson compuesto. Mini-ejemplo: en Simpson compuesto con $N$ paneles, el error se reduce en 16 veces al doblar $N$, nunca en 32.

5. Ignorar el error de redondeo por acumulación en mallas ultra-densas

Parece razonable asumir que usar $N = 10^8$ puntos en la regla del trapecio garantizará una precisión casi infinita. Falla porque sumar $10^8$ números de punto flotante en precisión simple (Float32) acumula un error de redondeo estocástico que supera con creces la ganancia del error de truncamiento, produciendo un resultado degradado por deriva numérica. Se detecta si el error calculado frente a una solución analítica deja de descender y comienza a crecer estocásticamente a partir de cierto $N$. Se corrige utilizando algoritmos de orden superior (Gauss-Legendre o Romberg) que alcanzan la precisión requerida con apenas decenas de puntos, o utilizando precisión Float64 con algoritmo de suma compensada de Kahan. Mini-ejemplo: sumar diez millones de franjas trapezoidales en Float32 acumula un error relativo de redondeo de $\sim 10^{-4}$, arruinando los 7 dígitos de mantisa teóricos.

Regla rápida de detección

• ¿$N$ es impar en Simpson? Inválido: añade un punto o usa la regla 3/8 en el último intervalo.
• ¿El integrando oscila rápido? Comprueba si $h < \pi / \omega$. Si no lo es, la integral es puro ruido espurio.

Aplicaciones físicas de la Integración Numérica (Cuadratura)

1. Electromagnetismo: cálculo del potencial electrostático de configuraciones continuas

Aparece en ingeniería de alta tensión y física de aceleradores al calcular el potencial electrostático $V(\mathbf{r})$ producido por electrodos metálicos cargados o haces de iones continuos en puntos arbitrarios del espacio tridimensional mediante la integral de Coulomb: $V(\mathbf{r}) = \frac{1}{4\pi\epsilon_0} \int_V \frac{\rho(\mathbf{r}')}{|\mathbf{r} - \mathbf{r}'|} d^3\mathbf{r}'$ Se calcula parametrizando la geometría del conductor y aplicando cuadraturas gaussianas bidimensionales o tridimensionales sobre elementos de contorno (Boundary Element Method, BEM). El resultado da la distribución espacial del potencial y permite derivar el campo eléctrico $\mathbf{E} = -\nabla V$ para predecir si se superará la rigidez dieléctrica del aire ($3 \text{ MV/m}$) provocando una descarga de arco eléctrico o efecto corona. El error peligroso es aproximar el potencial cuando el punto de observación $\mathbf{r}$ se ubica sobre la propia superficie cargada o muy cerca de ella; en ese límite el denominador $|\mathbf{r} - \mathbf{r}'| \to 0$ introduce una singularidad débil integrable que, si se evalúa con mallas regulares estándar, provoca picos numéricos artificiales de potencial infinito. Como control final se extrae la singularidad analíticamente mediante coordenadas polares locales (técnica de Duffy) y se comprueba la ley de Gauss en una superficie envolvente cerrada.

2. Termodinámica y astrofísica: potencia total radiada por el cuerpo negro (ley de Stefan-Boltzmann)

Aparece en astrofísica estelar para determinar la luminosidad bolométrica total de estrellas y galaxias integrando la radiancia espectral de Planck en todo el espectro electromagnético desde $\lambda = 0$ hasta $\lambda = \infty$: $B(T) = \int_0^\infty \frac{2 h c^2}{\lambda^5} \frac{1}{e^{\frac{h c}{\lambda k_B T}} - 1} d\lambda$ Se calcula transformando la integral a la variable adimensional $x = \frac{h c}{\lambda k_B T}$, lo que conduce a: $B(T) = \frac{2 k_B^4 T^4}{h^3 c^2} \int_0^\infty \frac{x^3}{e^x - 1} dx$ El resultado da la integral canónica $\int_0^\infty \frac{x^3}{e^x - 1} dx = \frac{\pi^4}{15}$, demostrando computacionalmente la ley de Stefan-Boltzmann $u = \sigma T^4$ con $\sigma = \frac{2\pi^5 k_B^4}{15 c^2 h^3} \approx 5.67037 \times 10^{-8} \text{ W/(m}^2\cdot\text{K}^4)$. El error peligroso es truncar el límite superior de integración en un valor bajo de $x$ (por ejemplo $x=5$) al evaluar filtros ópticos infrarrojos; cortar la integral prematuramente subestima la emisión térmica total en más de un $15\%$, sesgando el cálculo de la temperatura efectiva estelar. Como control final se verifica numéricamente mediante cuadratura de Gauss-Laguerre (diseñada para dominios semi-infinitos con peso $e^{-x}$) la convergencia al valor analítico $\pi^4/15 \approx 6.493939$.

3. Dinámica relativista y física de partículas: sección eficaz y factor de espacio fásico

Aparece en física de altas energías en colisionadores (CERN) al predecir la tasa de producción de partículas (como la producción del bosón de Higgs o la dispersión de Møller $e^- e^- \to e^- e^-$) integrando el elemento de matriz cuántico invariante de Feynman $|\mathcal{M}|^2$ sobre el espacio fásico de Lorentz de las partículas salientes: $\sigma = \frac{1}{2 s} \int |\mathcal{M}|^2 d\Phi_n$ Se calcula empleando algoritmos adaptativos de cuadratura multivariable como VEGAS (basado en Monte Carlo con refinamiento de malla por importancia). El resultado da la sección eficaz total $\sigma$ en picobarns ($1 \text{ pb} = 10^{-40} \text{ m}^2$), fijando cuántos eventos experimentales deben observarse por mes de operación del colisionador para confirmar un modelo de física más allá del Modelo Estándar. El error peligroso es aplicar cuadraturas estándar de paso fijo sobre resonancias cinemáticas de Breit-Wigner extremadamente estrechas ($|\mathcal{M}|^2 \propto \frac{1}{(s - M^2)^2 + M^2 \Gamma^2}$); si el ancho de desintegración $\Gamma$ es muy pequeño (como en el bosón $Z^0$, $\Gamma_Z \approx 2.5 \text{ GeV}$ frente a $M_Z \approx 91.2 \text{ GeV}$), una malla fija pasa por alto el pico de resonancia subestimando la producción de partículas en varios órdenes de magnitud. Como control final se realiza una transformación analítica de variable (técnica de multichannel flattening) que ensancha el pico de la resonancia antes de la integración numérica.

4. Física del estado sólido: densidad de estados y nivel de Fermi en metales

Aparece en física de la materia condensada al determinar el número total de electrones de conducción $n$ y la energía de Fermi $E_F(T)$ a temperatura finita integrando el producto de la densidad de estados $g(E)$ por la distribución de Fermi-Dirac: $n = \int_0^\infty g(E) f(E) dE = \int_0^\infty \frac{1}{2\pi^2} \left( \frac{2m^*}{\hbar^2} \right)^{3/2} \frac{\sqrt{E}}{e^{\frac{E - \mu}{k_B T}} + 1} dE$ Se calcula utilizando cuadratura de Gauss-Fermi o la expansión asintótica de Sommerfeld para bajas temperaturas ($k_B T \ll E_F$). El resultado da el potencial químico $\mu(T)$, la susceptibilidad magnética de Pauli y el calor específico electrónico lineal $C_e = \gamma T$. El error peligroso es evaluar la integral a temperaturas criogénicas ($T \to 0$) con paso fijo sin resolver la discontinuidad en escalón casi vertical que sufre la función de Fermi-Dirac en $E = \mu$; la derivada de $f(E)$ en ese punto es del orden de $1/(4 k_B T)$, lo que requiere mallas adaptativas con resolución sub-microelectrónvoltio alrededor del nivel de Fermi. Como control final se verifica que en el límite $T \to 0$ la integral converja exactamente al valor de Sommerfeld $n = \frac{2}{3} g(E_F) E_F$.

5. Mecánica de fluidos y aerodinámica: sustentación sobre perfiles alares por integración de presiones

Aparece en ingeniería aeronáutica y túneles de viento al calcular la fuerza total de sustentación $L$ y momento de cabeceo ejercido sobre un ala de avión integrando la distribución continua de presiones $P(s)$ y tensiones de cizalladura viscosa $\tau(s)$ a lo largo del contorno cerrado del perfil aerodinámico: $L = \oint_{\text{perfil}} [ -P(s) \cos\theta(s) - \tau(s) \sin\theta(s) ] ds$ Se calcula muestreando la presión mediante tomas estáticas piezométricas distribuidas en el intradós y trasdós del perfil y aplicando la regla del trapecio no uniforme compuesta sobre la curva curvilínea cerrada. El resultado da el coeficiente de sustentación $C_L$, certificando la velocidad de despegue y aterrizaje segura de una aeronave comercial. El error peligroso es utilizar una densidad de sensores insuficiente en el borde de ataque del ala; en ángulos de ataque elevados, el fluido se acelera violentamente alrededor del radio del borde de ataque generando un pico agudo de succión (depresión) extremadamente estrecho. Si no se densifica la integración en el primer $5\%$ de la cuerda alar, se pierde más del $30\%$ del área de sustentación real, falseando el punto de entrada en pérdida (stall). Como control final se comprueba mediante el teorema de Kutta-Joukowski ($L' = \rho_\infty V_\infty \Gamma$) que la sustentación integrada por presiones coincida rigurosamente con la circulación hidrodinámica $\Gamma = \oint \mathbf{v} \cdot d\mathbf{r}$ calculada en un contorno lejano.

Práctica guiada: Integración Numérica (Cuadratura)

Nivel 1: reconocimiento

Identifica la regla de cuadratura más adecuada para cada uno de los siguientes problemas físicos:

1. Calcular el flujo magnético $\Phi = \int_S \mathbf{B} \cdot d\mathbf{A}$ a partir de mediciones experimentales en 11 puntos equiespaciados tomadas por un sensor Hall sobre una placa plana.
2. Evaluar numéricamente la función de error de Gauss $\operatorname{erf}(x) = \frac{2}{\sqrt{\pi}} \int_0^x e^{-t^2} dt$ en una rutina de simulación cuántica con una precisión exigida de $10^{-14}$.
3. Integrar la energía potencial gravitatoria de un cuerpo heterogéneo donde los datos de densidad están tabulados en 8 puntos con espaciado irregular.
4. Integrar $\int_0^1 \frac{\cos(x)}{\sqrt{x}} dx$.

Solución:

1. Regla de Simpson 1/3 compuesta (11 puntos definen $N=10$ subintervalos, un número par perfecto para orden $\mathcal{O}(h^4)$).
2. Cuadratura de Gauss-Legendre de orden $n=8$ o $n=10$ (convergencia exponencial en funciones analíticas infinitamente diferenciables como $e^{-t^2}$).
3. Regla del Trapecio no uniforme compuesta (adapta automáticamente el área segmento a segmento mediante suma de trapecios).
4. Cambio de variable previo de regularización $x = u^2 \implies dx = 2u du$, transformando la integral en $2 \int_0^1 \cos(u^2) du$, la cual es completamente suave y libre de singularidades, resoluble por Gauss-Legendre o Simpson.

Nivel 2: aplicación directa

Aplica la regla de Simpson 1/3 compuesta con $N=2$ subintervalos ($h = 0.5$) para aproximar la integral de la energía cinética relativista acumulada: $I = \int_0^1 \frac{1}{1 + x^2} dx$ El valor analítico exacto es $\arctan(1) - \arctan(0) = \frac{\pi}{4} \approx 0.7853982$.

Solución:

1. Nodos: $a = 0, b = 1, N = 2 \implies h = 0.5$. $x_0 = 0.0, x_1 = 0.5, x_2 = 1.0$.
2. Valores del integrando $f(x) = \frac{1}{1 + x^2}$:
   • $f(x_0) = f(0) = \frac{1}{1 + 0} = 1.0000000$.
   • $f(x_1) = f(0.5) = \frac{1}{1 + (0.5)^2} = \frac{1}{1.25} = 0.8000000$.
   • $f(x_2) = f(1.0) = \frac{1}{1 + 1^2} = \frac{1}{2} = 0.5000000$.
3. Fórmula de Simpson 1/3: $I_{\text{Simp}} = \frac{h}{3} [f(x_0) + 4 f(x_1) + f(x_2)] = \frac{0.5}{3} [1.0000000 + 4(0.8000000) + 0.5000000]$ $= \frac{0.5}{3} [1.0 + 3.2 + 0.5] = \frac{0.5}{3} [4.7] = \frac{2.35}{3} \approx 0.7833333$
4. Comparación con el valor exacto: $E_{\text{abs}} = |0.7853982 - 0.7833333| = 0.0020649 \quad (\text{error de apenas el } 0.26\% \text{ con solo 3 puntos}).$

Nivel 3: elección de herramienta

Necesitas calcular la capacidad calorífica de un sólido cristalino según el modelo de Debye: $C_v(T) = 9 N k_B \left( \frac{T}{\Theta_D} \right)^3 \int_0^{\Theta_D/T} \frac{x^4 e^x}{(e^x - 1)^2} dx$ A temperaturas muy bajas ($T \to 0$), el límite superior $\Theta_D / T \to \infty$. ¿Cómo debes proceder numéricamente para evaluar la integral a bajas temperaturas sin perder precisión?

Solución: A temperaturas ultra-bajas, el límite superior tiende a infinito pero el integrando decae exponencialmente hacia cero como $x^4 e^{-x}$.

1. Analíticamente, la integral converge al límite infinito exacto: $\int_0^\infty \frac{x^4 e^x}{(e^x - 1)^2} dx = \frac{4\pi^4}{15} \approx 25.97576$
2. En lugar de forzar una cuadratura numérica con un límite grande arbitrario (como 1000) que introduciría cancelaciones y desbordamiento de $e^x$, se sustituye la integral por su valor analítico exacto $\frac{4\pi^4}{15}$, recuperando la clásica ley de Debye $T^3$: $C_v(T) \approx \frac{12\pi^4}{5} N k_B \left( \frac{T}{\Theta_D} \right)^3$
3. Para temperaturas intermedias ($0.1 < T/\Theta_D < 2$), se utiliza la cuadratura de Gauss-Legendre de 10 puntos en el intervalo acotado $[0, \Theta_D/T]$.

Nivel 4: interpretación física

Demuestra por qué la regla de Simpson 1/3, que se deduce integrando polinomios cuadráticos de segundo grado ($ax^2 + bx + c$), es de forma asombrosa exacta para cualquier polinomio cúbico de tercer grado ($f(x) = a_3 x^3 + a_2 x^2 + a_1 x + a_0$). ¿Qué propiedad de simetría física explica esta ganancia de orden gratuita?

Solución: Consideremos la integración de la función mónoma $f(x) = x^3$ en el intervalo simétrico $[-h, h]$:

1. La integral analítica exacta de una función impar en un intervalo simétrico es estrictamente cero: $\int_{-h}^h x^3 dx = \left[ \frac{x^4}{4} \right]_{-h}^h = \frac{h^4}{4} - \frac{h^4}{4} = 0$
2. Aplicamos la regla de Simpson con nodos en $-h, 0, h$: $S(h) = \frac{h}{3} [f(-h) + 4 f(0) + f(h)] = \frac{h}{3} [(-h)^3 + 4(0)^3 + h^3] = \frac{h}{3} [-h^3 + 0 + h^3] = 0$ El resultado de Simpson es exactamente cero, coincidiendo al $100\%$ con la integral exacta. Interpretación física: Por simetría de paridad, el error de truncamiento del término cúbico es una función impar respecto al punto medio del subintervalo. Al integrar sobre el intervalo completo, el error cometido en la mitad izquierda se cancela exactamente con el error de signo opuesto cometido en la mitad derecha. Como consecuencia, el término en $h^3$ desaparece del error y el primer término no nulo es el de orden cuatro $h^4$.

Nivel 5: error intencionado

Un estudiante calcula el potencial electrostático en el origen debido a una barra de carga uniforme de longitud $L=1$ con densidad lineal $\lambda$: $V = \frac{1}{4\pi\epsilon_0} \int_0^1 \frac{\lambda}{x} dx$ El estudiante programa la regla del trapecio con $N=100$ puntos: $x_0 = 0, x_1 = 0.01, \dots, x_{100} = 1.0$. El código aborta de inmediato con un error de ejecución por división entre cero. Explica el error conceptual físico y matemático del estudiante y cómo se regulariza el problema.

Fallo conceptual: El integrando $f(x) = 1/x$ posee una singularidad de polo simple en $x=0$.

1. Matemáticamente, la integral impropia diverge logarítmicamente: $\int_0^1 \frac{1}{x} dx = [\ln x]_0^1 = \ln(1) - \ln(0) = 0 - (-\infty) = +\infty$
2. Ningún método de cuadratura numérica puede integrar una singularidad divergente no integrable.
3. Físicamente, el origen $x=0$ coincide exactamente con el extremo de una barra de carga idealizada unidimensional con grosor cero. En la electrodinámica real, los conductores tienen un radio transversal finito $R > 0$. El potencial físico en el eje a una distancia de aproximación $R$ es: $V = \frac{\lambda}{4\pi\epsilon_0} \int_0^1 \frac{dx}{\sqrt{x^2 + R^2}} = \frac{\lambda}{4\pi\epsilon_0} \left[ \ln\left( x + \sqrt{x^2 + R^2} \right) \right]_0^1 = \frac{\lambda}{4\pi\epsilon_0} \ln\left( \frac{1 + \sqrt{1 + R^2}}{R} \right)$ Al incorporar el radio físico finito del conductor $R$, el integrando se vuelve estrictamente finito y suave en todo el intervalo, permitiendo su cálculo numérico exacto sin divisiones por cero.