Integración numérica

En análisis numérico, la integración​ numérica constituye una amplia gama de algoritmos para calcular el valor numérico de una integral definida y, por extensión, el término se usa a veces para describir algoritmos numéricos para resolver ecuaciones diferenciales. El término cuadratura numérica (a menudo abreviado a cuadratura) es más o menos sinónimo de integración numérica, especialmente si se aplica a integrales de una dimensión, a pesar de que para el caso de dos o más dimensiones (integral múltiple) también se utiliza.

Integración numérica empleando el método de Monte Carlo: Los nodos se distribuyen aleatoriamente en el intervalo de integración. Los puntos nuevos se representan en color azul marino, y los puntos antiguos en color azul claro. El valor de la integral tiende a 3,32

El problema básico considerado por la integración numérica es calcular una solución aproximada a la integral definida:

$\int _{a}^{b}f(x)\,dx$

Este problema también puede ser enunciado como un problema de valor inicial para una ecuación diferencial ordinaria, como sigue:

$y'(x)=f(x),\quad y(a)=0$

Encontrar y(b) es equivalente a calcular la integral. Los métodos desarrollados para ecuaciones diferenciales ordinarias, como el método de Runge-Kutta, pueden ser aplicados al problema reformulado. En este artículo se discuten métodos desarrollados específicamente para el problema formulado como una integral definida.

Razones para la integración numérica

Hay varias razones para llevar a cabo la integración numérica. La principal puede ser la imposibilidad de realizar la integración de forma analítica. Es decir, integrales que requerirían de un gran conocimiento y manejo de matemática avanzada pueden ser resueltas de una manera más sencilla mediante métodos numéricos. Incluso existen funciones integrables pero cuya primitiva no puede ser calculada, siendo la integración numérica de vital importancia. La solución analítica de una integral nos arrojaría una solución exacta, mientras que la solución numérica nos daría una solución aproximada. El error de la aproximación, que depende del método que se utilice y de qué tan preciso sea, puede llegar a ser tan pequeño que es posible obtener un resultado idéntico a la solución analítica en las primeras cifras decimales.

Métodos para integrales unidimensionales

Los métodos de integración numérica pueden ser descritos generalmente como combinación de evaluaciones del integrando para obtener una aproximación a la integral. Una parte importante del análisis de cualquier método de integración numérica es estudiar el comportamiento del error de aproximación como una función del número de evaluaciones del integrando. Un método que produce un pequeño error para un pequeño número de evaluaciones es normalmente considerado superior. Reduciendo el número de evaluaciones del integrando se reduce el número de operaciones aritméticas involucradas, y por tanto se reduce el error de redondeo total. También, cada evaluación cuesta tiempo, y el integrando puede ser arbitrariamente complicado.

Cabe decir que el problema sólo tiene interés si no es posible conocer con exactitud el valor del error. El error por definición es la diferencia entre el valor aproximado y el valor exacto. En los casos en los que se puede conocer el error, basta con añadirlo al valor aproximado para obtener el valor exacto. Pero en estos casos no es necesaria una integración numérica. Por tanto los métodos trabajan o bien con metas superiores del error o bien con valores estimados del error.

Entre dos métodos que tienen el mismo error, se considera superior el que necesita un menor número de evaluaciones de la función.

El hecho de reducir el número de evaluaciones de la función integrando reduce el número de operaciones aritméticas implicadas. Esto reduce el tiempo de cálculo. En los casos en que los cálculos se realizan con operaciones de coma flotante (que son la mayoría) el hecho de reducir el número de operaciones también reduce el error de redondeo total.

De todos modos, un modo de integración por «fuerza bruta» puede hacerse siempre, de un modo muy simplista, evaluando el integrando con incrementos muy pequeños.

Métodos basados en funciones de interpolación

Hay una extensa familia de métodos que se basan en aproximar la función a integrar $f(x)$ por otra función $g(x)$ de la cual se conoce la integral exacta. La función que sustituye la original se encuentra de forma que en un cierto número de puntos tenga el mismo valor que la original. Como los puntos extremos forman parte siempre de este conjunto de puntos, la nueva función se llama una interpolación de la función original. Cuando los puntos extremos no se utilizan para encontrar la función que sustituye a la original entonces se dice extrapolación. Típicamente estas funciones son polinomios.

Fórmulas de Newton-Cotes

La interpolación con polinomios evaluada en puntos igualmente separados en $\ [a,b]$ da las fórmulas de Newton-Cotes, de las que la regla del rectángulo, la del trapecio y la de Simpson son ejemplos. Si se escogen los nodos hasta $k=n+1$ incluidos los extremos del rango será la fórmula de Newton-Cotes cerrada y si se escogen $k=n-1$ o que no incluyen en los datos tabulados los extremos será la fórmula de Newton-Cotes abierta.

Regla del rectángulo

El método más simple de este tipo es hacer a la función interpoladora ser una función constante (un polinomio de orden cero) que pasa a través del punto $(a,f(a))$. Este método se llama la regla del rectángulo:

$\int _{a}^{b}f(x)dx\sim (b-a)f(a)$

Regla del punto medio

Ilustración de la regla del punto medio.

Si en el método anterior la función pasa a través del punto $\left({\frac {a+b}{2}},f\left({\frac {a+b}{2}}\right)\right)$ este método se llama la regla del punto medio:

$\int _{a}^{b}f(x)dx\sim (b-a)\,f\!\left({\frac {a+b}{2}}\right)$

Regla del trapecio

Ilustración de la regla del trapecio

Corresponde al caso donde el polinomio sustituto de la ecuación original es del primer orden aunque la fórmula original no pase por el punto $\left({\frac {a+b}{2}},f\left({\frac {a+b}{2}}\right)\right)$. Entonces si:

$I=\int _{a}^{b}\left[f(a)+{\frac {f(b)-f(a)}{b-a}}(x-a)\right]dx$

la regla de integración es:

$I=(b-a){\frac {f(a)+f(b)}{2}}$

que viene a ser la misma fórmula que la usada en la regla del punto medio. ​

Regla de Simpson

Ilustración de la regla de Simpson.

La función interpoladora puede ser un polinomio de grado 2 que pasa a través de los puntos $\ (a,f(a))$, $({\frac {a+b}{2}},f({\frac {a+b}{2}}))$ y $\ (b,f(b))$. Este método se llama regla de Simpson:

$\int _{a}^{b}f(x)\,dx\approx {\frac {b-a}{6}}\left[f(a)+4f\left({\frac {a+b}{2}}\right)+f(b)\right]$.

Reglas compuestas

Para cualquier regla interpoladora, se puede hacer una aproximación más precisa dividiendo el intervalo $\ [a,b]$ en algún número $\ n$ de subintervalos, hallando una aproximación para cada subintervalo, y finalmente sumando todos los resultados. Las reglas que surgen de hacer esto se llaman reglas compuestas, y se caracterizan por perder un orden de precisión global frente a las correspondientes simples, si bien globalmente dan valores más precisos de la integral, a costa eso sí de incrementar significativamente el coste operativo del método. Por ejemplo, la regla del trapecio compuesta puede expresarse como:

$\int _{a}^{b}f(x)dx\sim {\frac {b-a}{n}}\left({\frac {f(a)+f(b)}{2}}+\sum _{k=1}^{n-1}f\left(a+k{\frac {b-a}{n}}\right)\right)$

donde los subintervalos tienen la forma $\ [kh,(k+1)h]$ con
$\ h={\frac {b-a}{n}}$
y $\ k=0,1,2,\ldots ,n-1$.

Métodos de extrapolación

La precisión de un método de integración del tipo Newton-Cotes es generalmente una función del número de puntos de evaluación. El resultado es usualmente más preciso cuando el número de puntos de evaluación aumenta, o, equivalentemente, cuando la anchura del paso entre puntos decrece. ¿Qué pasa cuando la anchura del paso tiende a cero? Esto puede responderse extrapolando el resultado de dos o más anchuras de paso (extrapolación de Richardson). La función de extrapolación puede ser un polinomio o una función racional. Los métodos de extrapolación están descritos en más detalle por Stoer y Bulirsch (Sección 3.4). En particular, al aplicar el método de extrapolación de Richardson a la regla del trapecio compuesta se obtiene el método de Romberg.

Cuadratura de Gauss

Si se permite variar los intervalos entre los puntos de interpolación, se encuentra otro grupo de fórmulas de integración, llamadas fórmulas de cuadratura de Gauss. Una regla de cuadratura de Gauss es típicamente más precisa que una regla de Newton-Cotes que requiera el mismo número de evaluaciones del integrando, si el integrando es suave (es decir, si se puede derivar muchas veces).

Algoritmos adaptativos

Si f no tiene muchas derivadas definidas en todos sus puntos, o si las derivadas toman valores muy elevados, la integración gausiana es a menudo insuficiente. En este caso, un algoritmo similar al siguiente lo haría mejor:

def integral(f, a, b):
    """Este algoritmo calcula la integral definida de una función
    en el intervalo [a,b], adaptativamente, eligiendo pasos más
    pequeños cerca de los puntos problemáticos.
    h0 es el paso inicial."""
    
    x = a
    h = h0
    acumulador = 0
    while x < b:
        if x+h > b:  h = b - x
        if error de la cuadratura sobre [x,x+h] para f es demasiado grande:
            haz h más pequeño
        else:
            acumulador += integral(f, x, x+h)
            x += h
            if error de la cuadratura sobre [x,x+h] es demasiado pequeño:
                haz h más grande
    return acumulador

Algunos detalles del algoritmo requieren mirarlo con cuidado. Para muchos casos, estimar el error de la integral sobre un intervalo para un función f no es obvio. Una solución popular es usar dos reglas de integración distintas, y tomar su diferencia como una estimación del error de la integral. El otro problema consiste en decidir qué es «demasiado grande» o «demasiado pequeño». Un criterio posible para «demasiado grande» es que el error de la integral no sea mayor que th, donde t, un número real, es la tolerancia que queremos tener para el error global. Pero también, si h es ya minúsculo, puede no valer la pena hacerlo todavía más pequeño si el error de la integral es aparentemente grande. Este tipo de análisis de error usualmente se llama «a posteriori» ya que calculamos el error después de haber calculado la aproximación.

La heurística para integración adaptativa está discutida en Forsythe et al (sección 5.4).

Estimación del error conservativa (a priori)

Supongamos que $\ f$ tiene una primera derivada sobre $\ [a,b]$ acotada. El teorema del valor medio para $\ f$, para $\ x<b$, da

$\ (x-a)f'(y_{x})=f(x)-f(a)$

para algún $\ y_{x}$ en $\ [a,x]$ dependiendo de $\ x$. Si integramos en $\ x$ de $\ a$ a $\ b$ en ambos lados de la igualdad y tomamos valores absolutos, tenemos

$\left|\int _{a}^{b}f(x)dx-(b-a)f(a)\right|=\left|\int _{a}^{b}(x-a)f'(y_{x})dx\right|$

Se puede aproximar más la integral en el lado derecho metiendo el valor absoluto en el integrando, y reemplazando el término en $\ f$ por una cota superior:

$\left|\int _{a}^{b}f(x)dx-(b-a)f(a)\right|\leq {\frac {(b-a)^{2}}{2}}\sup _{a\leq x\leq b}\left|f'(x)\right|$

Así, si aproximamos la integral $\int _{a}^{b}f(x)dx$ por su regla de integración $\ (b-a)f(a)$, el error no es mayor que el lado derecho de la ecuación.

Integrales múltiples

Los métodos de integración que se han comentado hasta aquí se han diseñado todos para calcular integrales de una dimensión.

Para calcular integrales de diversas dimensiones, un enfoque es expresar la integral múltiple como repetición de integrales de una dimensión haciendo uso del teorema de Fubini.

Este enfoque lleva a una cantidad de evaluaciones de la función que crece exponencialmente a medida que crece el número de dimensiones. Se conocen dos métodos para superar esta llamada maldición de la dimensión.

Montecarlo

Los métodos de Montecarlo y métodos de cuasi-Montecarlo son fáciles de aplicar a integrales multidimensionales, y pueden producir una mejor exactitud por el mismo número de evaluaciones de la función que en integraciones repetidas empleando métodos unidimensionales. Una clase grande de métodos útiles de Montecarlo son los llamados algoritmos de Cadena de Markov de Montecarlo, los cuales incluyen el algoritmo de Metropolis-Hastings y muestreo de Gibbs.

Parrillas dispersas

Las parrillas dispersas fueron desarrolladas inicialmente por Smolyak para la cuadratura de funciones de muchas dimensiones. El método se basa siempre en un método de cuadratura unidimensional, pero realiza una más sofisticada combinación de resultados de cada variable.

Conexión con las ecuaciones diferenciales

El problema de evaluar la integral

$\int _{a}^{b}f(x)\,dx$

se puede reducir a un problema de valor inicial para una ecuación diferencial ordinaria. Si la integral de más arriba se escribe como I(b), entonces la función I satisface

$I'(x)=f(x),\quad I(a)=0.$

Los métodos desarrollados para resolver ecuaciones diferenciales ordinarias, como por ejemplo los métodos de Runge-Kutta, se pueden aplicar para resolver el problema y por tanto se pueden utilizar para evaluar la integral.

Programas para integración numérica

La integración numérica es uno de los problemas estudiados más intensivamente en el análisis numérico. Entre las muchas implementaciones en programas se encuentran:

• QUADPACK (parte de SLATEC) (código fuente): QUADPACK es una colección de algoritmos en Fortran para integración numérica basada en reglas gausianas.
• GSL: GNU Scientific Library. La biblioteca Científica de GNU (GSL) es una biblioteca numérica escrita en C que provee una amplia gama de rutinas matemáticas, como la integración por Montecarlo.
• ALGLIB: Es una colección de algoritmos en C# / C++ / Delphi / Visual Basic / etc., para la integración numérica.

Se pueden encontrar algoritmos de integración numérica en GAMS class H2.

Bibliografía

• Philip J. Davis and Philip Rabinowitz, Methods of Numerical Integration.
• George E. Forsythe, Michael A. Malcolm, and Cleve B. Moler. Computer Methods for Mathematical Computations. Englewood Cliffs, NJ: Prentice-Hall, 1977. (Ver Capítulo 5.)
• William H. Press, Brian P. Flannery, Saul A. Teukolsky, William T. Vetterling. Numerical Recipes in C. Cambridge, UK: Cambridge University Press, 1988. (Ver Capítulo 4.)
• Josef Stoer and Roland Bulirsch. Introduction to Numerical Analysis. New York: Springer-Verlag, 1980. (Ver Capítulo 3.)
• Jon M. Smith Recent Developments in Numerical Integration, J. Dynam. Sys., Measurement and Control 96, Ser. G-1, No. 1, 61-70, Mar. 1974.

Enlaces externos

• Wikimedia Commons alberga una categoría multimedia sobre Integración numérica.
• Integración: Antecedentes, Simulaciones, etc. en el Holistic Numerical Methods Institute

• Datos: Q753445
• Multimedia: Numerical integration / Q753445