Cuadratura explicada: métodos numéricos para aproximar integrales

Hay una categoría de integrales definidas que se resisten al cálculo a mano. No porque lo estés haciendo mal. Algunas funciones, como e^(−x²) o sin(x)/x, no tienen primitivas en forma cerrada. Cuando pasa eso, necesitas otro enfoque: evaluar la función en unos cuantos puntos bien elegidos y combinar esos valores para estimar el área bajo la curva. A eso lo llamamos cuadratura.

Este artículo explica qué es la cuadratura, cómo funciona el marco general y los métodos principales, desde la intuitiva regla del punto medio hasta el sorprendentemente preciso enfoque gaussiano. También veremos métodos compuestos, fuentes de error, estrategias adaptativas y cómo elegir el método adecuado para el problema que tienes delante.

¿Qué es la cuadratura?

La palabra "cuadratura" es antigua. Viene del problema geométrico clásico de construir un cuadrado con la misma área que una figura dada, que no es más que otra forma de decir "calcula el área". Hoy el término se usa para la aproximación numérica de una integral definida.

La idea central es un cambio de enfoque: en lugar de hallar una primitiva exacta, muestreas la función f(x) en un conjunto de puntos llamados nodos, multiplicas cada valor por un peso asociado y sumas los resultados:

Los términos cuadratura e integración numérica significan esencialmente lo mismo. Cuadratura se usa más en contextos matemáticos; integración numérica es el término habitual en computación científica.

Cómo funciona la cuadratura numérica

Cada regla de cuadratura opera en un intervalo [a, b]. Dentro de ese intervalo, selecciona un conjunto de nodos xᵢ y asigna a cada uno un peso wᵢ. La aproximación es la suma ponderada de arriba. Lo que cambia entre métodos es dónde van los nodos y cómo se asignan los pesos. Nada más. La regla del trapecio, la regla de Simpson y la cuadratura gaussiana son variaciones sobre este mismo patrón.

El error de aproximación es la diferencia entre la suma ponderada y la integral verdadera. Todos los métodos introducen algo de error, y lo rápido que se reduce al añadir más nodos (o acotar más los subintervalos) es lo que separa un buen método de uno excelente. Quédate con esta idea: elegir un método de cuadratura es, en realidad, decidir cómo gastar tus evaluaciones de la función.

Métodos de cuadratura comunes

Los métodos siguientes van de lo más intuitivo a lo más sofisticado.

Comparativa visual de las reglas de integración de punto medio, trapecio y Simpson. Imagen del autor.

Regla del punto medio

Toma el punto medio de [a, b], evalúa f ahí y multiplícalo por el ancho (b − a). La suposición es que el valor de la función en el centro representa el promedio en todo el intervalo.

Suena casi demasiado simple, y en un único intervalo lo es. Pero aquí viene la primera sorpresa: la regla del punto medio es exacta para funciones constantes y lineales, y su error escala como O(h³) por intervalo, donde h = b − a. De hecho, eso es mejor que la regla del trapecio en funciones suaves, algo que suele pillar desprevenido a más de uno.

Regla del trapecio

Usa ambos extremos en lugar del centro. Conéctalos con una recta y calcula el área del trapecio resultante: (b − a) · [f(a) + f(b)] / 2. Es exacta para funciones lineales, con error O(h³) por intervalo, del mismo orden que la regla del punto medio, aunque con una constante mayor.

Probablemente sea el método de cuadratura más enseñado, lo que le ha dado una fama desmesurada. Para funciones bien comportadas es razonable. Pero la versión compuesta de la regla de Simpson casi siempre rinde mejor con el mismo número de evaluaciones.

Regla de Simpson

Ajusta un polinomio cuadrático por tres puntos: los dos extremos y el punto medio. Integra ese polinomio exactamente:

La evaluación extra en el centro cuesta una llamada adicional a la función, pero compensa: la regla de Simpson es exacta para polinomios de hasta grado tres, y el error por intervalo escala como O(h⁵), dos órdenes mejor que la del trapecio. Para funciones suaves, suele ser el primer método al que merece la pena recurrir.

Cuadratura gaussiana

Aquí es donde la cosa se pone realmente interesante. Las reglas de punto medio, trapecio y Simpson usan nodos equiespaciados. La cuadratura gaussiana plantea otra pregunta: ¿y si eligiéramos tanto la ubicación de los nodos como los pesos para maximizar la precisión?

Con n nodos, la cuadratura gaussiana (concretamente la variante de Gauss-Legendre) puede integrar exactamente cualquier polinomio de grado hasta 2n − 1. Dos nodos te dan exactitud para cúbicos; cinco nodos cubren polinomios de grado 9. Los nodos son las raíces de los polinomios de Legendre y se agrupan hacia los extremos del intervalo en lugar de ser equiespaciados. La precisión extra viene de permitir que los nodos se muevan, no de dónde caen por casualidad. Su concentración cerca de los extremos es un subproducto de las matemáticas, y más abajo veremos por qué.

Así es como se comparan los cuatro métodos:

Métodos de cuadratura compuesta

Aplicar una regla una sola vez en todo el intervalo [a, b] suele no ser suficiente. La solución estándar: divide el intervalo en n subintervalos iguales de anchura h = (b − a)/n y aplica la regla básica en cada uno.

Reglas compuestas de punto medio y trapecio

Para la regla del trapecio compuesta, evalúa f en n+1 puntos equiespaciados y suma n trapecios. El error total baja a O(h²). Si reduces h a la mitad, el error cae por un factor de cuatro. La regla compuesta del punto medio sigue la misma lógica y tiene el mismo orden de error global.

Regla de Simpson compuesta

La versión compuesta requiere un número par de subintervalos. Evalúa f en n+1 puntos y obtiene un error global O(h⁴). Al reducir h a la mitad, el error baja por un factor de dieciséis. Por eso la regla de Simpson compuesta suele ser la opción por defecto cuando necesitas una precisión moderada en una función suave.

El siguiente código en Python muestra cómo mejora la aproximación al añadir subintervalos para ∫[0 a π] sin(x) dx (valor exacto: 2.0):

import numpy as np

def composite_simpson(func, lower, upper, num_intervals):
    """Approximate integral using composite Simpson's rule."""
    if num_intervals % 2 != 0:
        raise ValueError("num_intervals must be even")
    nodes = np.linspace(lower, upper, num_intervals + 1)
    values = func(nodes)
    h = (upper - lower) / num_intervals
    weights = np.ones(num_intervals + 1)
    weights[1:-1:2] = 4   # odd-indexed interior points
    weights[2:-2:2] = 2   # even-indexed interior points
    return (h / 3) * np.dot(weights, values)

exact = 2.0
for n in [4, 8, 16, 32]:
    approx = composite_simpson(np.sin, 0, np.pi, n)
    print(f"n={n:3d}: approx={approx:.10f}, error={abs(approx - exact):.2e}")

Con solo 32 intervalos, el error baja de 10⁻⁶. Eso es convergencia O(h⁴).

Cuadratura gaussiana

La idea detrás de la cuadratura gaussiana merece más que una fila en una tabla.

Por qué los nodos equiespaciados no siempre son óptimos

Cuando los nodos están equiespaciados, estás gastando tu presupuesto de evaluaciones de una manera que puede no ajustarse a la dificultad del integrando. Hay regiones que necesitan más resolución; en otras se desperdician evaluaciones en tramos planos. La cuadratura gaussiana evita esto eligiendo nodos que son demostrablemente óptimos para el grado de polinomio que puedes integrar exactamente.

Para n nodos, la regla de Gauss-Legendre selecciona como nodos las n raíces del n-ésimo polinomio de Legendre Pₙ(x) y luego calcula los pesos para que cualquier polinomio de grado ≤ 2n − 1 se integre sin error. En la práctica, numpy.polynomial.legendre.leggauss(n) devuelve ambos. Los nodos se agrupan hacia los extremos de [−1, 1], el intervalo de referencia estándar, donde el error de interpolación polinómica tiende a ser mayor y donde compensa aumentar la densidad de muestreo.

Cuándo usarla

La cuadratura gaussiana brilla cuando tienes una función suave y analítica y necesitas alta precisión con pocas evaluaciones. Cinco puntos de Gauss-Legendre pueden superar a la regla de Simpson compuesta con decenas de intervalos en integrandos bien comportados. Pero no es adecuada para funciones con esquinas, discontinuidades o datos disponibles solo en puntos equiespaciados fijos. En esos casos, la cuidadosa colocación de nodos juega en tu contra.

Error y precisión en cuadratura

Ningún método numérico se libra del error. Dos fuentes distintas afectan a todo resultado de cuadratura, y entender ambas es clave para saber cuándo fiarte de tu respuesta.

Las fuentes de error incluyen:

• Error de discretización (también llamado error de truncamiento) es la brecha entre la integral exacta y lo que devuelve tu regla de cuadratura en aritmética perfecta. Esto es lo que capturan las fórmulas O(hᵏ). Se reduce al añadir más puntos o disminuir la anchura de los intervalos.
• Error de coma flotante es un problema distinto. Sumar muchas cantidades pequeñas en aritmética de precisión finita introduce redondeos que no disminuyen con h y pueden incluso crecer. Con reglas compuestas usando miles de subintervalos, puedes llegar a un suelo donde domina el ruido numérico y añadir intervalos deja de ayudar. Para la mayoría de tareas de computación científica, alcanzarás tu tolerancia deseada mucho antes de ese suelo. Pero existe.

La suavidad de la función, el ancho del intervalo y la elección del método influyen en el error de discretización. Una función con rasgos bruscos necesita más intervalos o un enfoque adaptativo.

Cuadratura adaptativa

Una malla uniforme desperdicia evaluaciones cuando el integrando es suave en casi todo el dominio pero complicado en unos pocos puntos.

Cómo funciona

Aplica una regla de cuadratura a un intervalo, luego aplica una regla de orden superior (o la misma regla en las dos mitades) y compara. La diferencia estima el error local. Si está por debajo de la tolerancia, acepta el intervalo. Si no, divídelo y repite recursivamente en cada mitad. Repite hasta que cada subintervalo cumpla la tolerancia.

La ganancia de eficiencia sobre el refinamiento uniforme puede ser grande cuando el integrando tiene rasgos aislados: un pico estrecho, un quiebro, una oscilación rápida rodeada de zonas suaves. Los métodos adaptativos concentran las evaluaciones justo donde la función las necesita.

scipy.integrate.quad()

En Python, scipy.integrate.quad() implementa cuadratura gaussiana adaptativa (concretamente el algoritmo DQAGSE de QUADPACK). Pasa tu función y los límites; devuelve la estimación y una cota del error. Para la mayoría de las integrales unidimensionales, es la herramienta adecuada, sin necesidad de ajustar nada a mano.

from scipy import integrate
import numpy as np

# Integrate a function with a peak near x=0.5
def peaked_func(x):
    return np.exp(-100 * (x - 0.5)**2)

result, error_bound = integrate.quad(peaked_func, 0, 1)
print(f"Estimate: {result:.10f}")
print(f"Error bound: {error_bound:.2e}")

quad() resuelve este pico estrecho sin problema. Una Simpson uniforme necesitaría cientos de subintervalos para igualarlo.

Cómo elegir un método de cuadratura

Empieza con scipy.integrate.quad() para integrales unidimensionales en Python salvo que tengas un motivo para no hacerlo. Es adaptativo, preciso y muy probado. Todo lo que sigue es relevante cuando implementas algo tú mismo, trabajas en un lenguaje sin una buena librería adaptativa o estás optimizando para un caso concreto.

• Funciones suaves y analíticas con fórmula conocida: la cuadratura gaussiana ofrece la mayor precisión por evaluación. Si buscas precisión moderada y simplicidad, la regla de Simpson compuesta es una apuesta segura.
• Funciones con quiebros o discontinuidades: los métodos equiespaciados (trapecio, Simpson) tienden a superar a la cuadratura gaussiana aquí, porque la colocación fija de nodos no se ve afectada por la irregularidad. Divide el intervalo de integración en la discontinuidad conocida y trata cada tramo por separado.
• Datos en puntos fijos y equiespaciados: la cuadratura gaussiana no es opción. No puedes evaluar la función en nodos arbitrarios, así que usa trapecio compuesto o Simpson compuesto.
• Integrales de mayor dimensión: la cuadratura sufre pronto la maldición de la dimensionalidad. Con n puntos por dimensión y d dimensiones, necesitas nᵈ evaluaciones. Diez dimensiones con diez puntos cada una son 10¹⁰ evaluaciones, inviable. Monte Carlo escala mucho mejor por encima de unas cinco dimensiones, así que ahí es donde conviene cambiar.

Cuadratura frente a métodos numéricos relacionados

La cuadratura no existe en el vacío. Ayuda ver dónde encaja respecto a otras herramientas que probablemente estés valorando.

Cuadratura vs. integración analítica

La integración analítica obtiene una primitiva exacta; la cuadratura aproxima un número. Si existe una primitiva exacta y puedes encontrarla, úsala. La cuadratura se gana su sitio cuando la integración exacta es imposible o poco práctica.

Cuadratura vs. integración Monte Carlo

Ambas aproximan integrales definidas, pero de forma distinta. La cuadratura es determinista: evalúa f en nodos fijos y preelegidos, con convergencia ligada a la suavidad de la función y al número de nodos. Monte Carlo toma muestras aleatorias y las promedia, con error que escala como O(1/√N) independientemente de la dimensión. Es una convergencia más lenta en baja dimensión, pero no depende de la dimensión, lo que la convierte en la opción práctica a partir de unas cinco dimensiones.

Cuadratura vs. interpolación

Algunas reglas de cuadratura, incluida la de Simpson, funcionan ajustando un polinomio a los puntos de evaluación e integrándolo exactamente. Pero el objetivo difiere de la interpolación. La interpolación reconstruye la función; la cuadratura solo quiere un número. El polinomio es una herramienta, no el resultado.

Conclusión

Lo satisfactorio de la cuadratura es cuánto se aprovecha una única idea organizadora. Cada método aquí, punto medio, trapecio, Simpson, gaussiano, no es más que una respuesta distinta a la misma pregunta: dado un presupuesto de evaluaciones de la función, ¿dónde las pones y cómo las ponderas?

Ajustar mejores polinomios a los puntos de evaluación reduce el error rápidamente. Pasar a nodos elegidos de forma óptima exprime aún más precisión con el mismo presupuesto. Los métodos adaptativos rematan el conjunto, enviando las evaluaciones justo a donde la función realmente las necesita.

En Python, scipy.integrate.quad() es la mejor opción por defecto para problemas 1D. Todo lo demás en este artículo sirve para entender qué pasa por debajo del capó o cuando vas a implementarlo tú mismo. Si quieres seguir desarrollando habilidades numéricas con SciPy, haz nuestro curso Introduction to Optimization in Python.