PDF de programación - Integración adaptativa con quadpack

Integración adaptativa con quadpack

Pablo Santamaría
v0.1 (Octubre 2011)

Planteo del problema

Supongamos que se desea calcular una estimación numérica Q de la integral definida

Z b

a

I =

f(x) dx

tal que el error |E| = |Q − I| sea menor o igual que cierta tolerancia
prefijada. Si la función f se
puede evaluar en cualquier punto del intervalo de integración podríamos comenzar considerando
alguna fórmula de cuadratura numérica compuesta como ser la regla compuesta del trapecio sobre
n subintervalos igualmente espaciados una distancia h = (b − a)/n,

f(x0) + f(xn) +

Q = h
2

f(xi),

n−1X

i=1

siendo x0 = a, xn = b, xi = x0 + ih para i = 0, 1, . . . , n, cuyo término de error es

para algún ξ ∈ (a, b). Si disponemos de una cota M de f00 entonces la desigualdad

E = − h2(b − a)

12

f00(ξ),

|E| ≤ h2(b − a)

12 M ≤
,

permite determinar el paso h (y el valor de n) de una aproximación Q que efectivamente está
dentro de la tolerancia prefijada. Sin embargo, para una implementación algorítmica, debemos
encontrar una manera de estimar la precisión de la aproximación Q que no involucre determinar
f00. Para hacer esto utilizamos un principio básico del análisis numérico: la estimación del error de
una aproximación resulta de comparar con una aproximación más precisa. Supongamos, entonces,
que además de Q disponemos de otra aproximación ˆQ de I que fue obtenida con un método más
preciso. La ecuación

E = I − Q = ˆQ − Q + (I − ˆQ),

nos dice que cuando ˆQ es más precisa que Q, el error en Q puede estimarse como ˆQ − Q. Ahora
bien, el cálculo de ˆQ agrega un costo computacional que debería llevarse al mínimo reutilizando
para su cómputo las evaluaciones de f ya utilizadas para el cálculo de Q. Por ejemplo, en nuestro
caso, donde utilizamos la regla del trapecio para evaluar Q podemos utilizar la regla compuesta
de Simpson para determinar ˆQ (y por tanto, estimar el error de Q), dado que sabemos que dicho
método es considerablemente más preciso que el primero y además podemos utilizar los mismos
puntos xi (en tanto n + 1 sea un número impar).

Disponiendo de una fórmula de cuadratura compuesta y la posibilidad, según lo anterior,
de estimar el error cometido, se puede entonces implementar un procedimiento automático de

1

integración numérica como sigue: Comenzando con una paso de integración h (digamos h = b− a)
calculamos la aproximación Q y la estimación de su error. Si dicha estimación del error es menor
o igual que la tolerancia
prefijada entonces Q es la aproximación buscada a I. De lo contrario,
dividimos el paso h de integración a la mitad y comenzamos de nuevo. Aunque funciona, este
procedimiento es, en general, ineficiente debido a la gran cantidad de evaluaciones del integrando
que deben realizarse al utilizar un paso de integración constante sobre todo el intervalo [a, b] en
cada iteración. Un mejor esquema de integración surge teniendo presente que en aquellas regiones
del intervalo de integración donde el integrando f no presenta grandes variaciones se pueden
tomar subintervalos “grandes” y que en regiones donde la variación de f es “grande” podemos
tomar subintervalos “pequeños” de manera tal que la suma de las contribuciones al error total sea
menor que la tolerancia prefijada. Los procedimientos de integración que adaptan la longitud de los
subintervalos al comportamiento del integrando se llaman, apropiadamente, métodos adaptativos
de cuadratura.

Un método adaptativo de cuadratura procede como sigue: se aplica la regla de cuadratura sobre
todo el intervalo (esto es, con h = b−a) para obtener una aproximación Q a la integral y se estima su
error (a través de la aplicación de una regla de cuadratura de mayor precisión). Si el error estimado
es mayor que la tolerancia prefijada, se divide el intervalo de integración a la mitad y se aplica la
regla de cuadratura a cada uno de los subintervalos y se estima sus respectivos errores. Ahora, si la
suma de dichas estimaciones de los errores excede la tolerancia requerida, entonces el subintervalo
con mayor cota de error es, a su vez, subdividido a la mitad y así siguiendo hasta que la estimación
del error total sea menor que la tolerancia prefijada. Con más precisión, en una dada etapa del
procedimiento el intervalo [a, b] se encuentra particionado en subintervalos dados por los puntos a =
α1 < β1 = α2 < β2 = α3 < ··· < βn = b donde se tiene una estimación Qi para cada integral Ii de
f sobre [αi, βi] y una estimación Ei del error de Qi. A partir de estos valores se tiene una estimación

Q =P Qi para la integral I con una estimación E =P Ei de su error. Si esta aproximación Q
no está dentro de la tolerancia prefijada, esto es, si no se cumple la condición |E| = |P Ei| ≤
,
sus cantidades asociadas. El procedimiento es repetido hasta que |E| = |P Ei| ≤
1. Al final

entonces el subintervalo [αi, βi] con mayor cota |Ei| de error es escogido para subdividirlo. Dicho
intervalo es entonces particionado en dos mitades: [αi, (αi+βi)/2], [(αi+βi)/2, βi] y se calculan las
aproximaciones a las integrales sobre cada uno de tales subintervalos y sus respectivas estimaciones
de error. Estos dos subintervalos y las cantidades asociadas reemplazan al subintervalo [αi, βi] y

del proceso, la estimación Q de la integral ha sido calculada, dentro de la tolerancia prefijada,
aplicando la regla de cuadratura en una partición del intervalo [a, b] dada por subintervalos [αi, βi]
que tienden a ser más pequeños en aquellas regiones donde f presenta mayores variaciones.

Al igual que en nuestra primera aproximación naïve, un método adaptivo requiere de una
fórmula de cuadratura y una estimación de el error cometido a través de una regla de cuadratura
de mayor precisión que, con el fin de mantener mínimo el esfuerzo computacional, debería reutilizar
los mismos puntos que la regla inicial. La combinación de la regla de trapecios y Simpson parecen
buenos candidatos y, con más generalidad, se pueden utilizar reglas de cuadratura de Newton-
Cotes de diferentes órdenes, puesto que, al tener un paso igualmente espaciado compartirán puntos
en común. Sin embargo, sabemos que las reglas de cuadratura gaussianas son más precisas que
las reglas de Newton-Cotes para el mismo número de puntos. Pero, desafortunadamente, reglas
gaussianas de diferente órden no tienen puntos en común, con lo cual, la estimación del error con
reglas gaussianas requiere la evaluación del integrando f sobre todos los conjuntos de puntos de
ambas reglas. Para evitar tal trabajo computacional se han desarrollado las reglas de cuadratura de
Gauss-Kronrod. Una dada regla de Gaus-Kronrod consiste en una regla de Gauss Gn de n-puntos
y una regla de Kronrod K2n+1 de (2n+1)-puntos que han sido escogidos de manera óptima sujetos
a la condición de que todos los puntos de Gn son reusados en K2n+1

con lo cualP|Ei| ≤
. Dado que |P Ei| ≤P|Ei| el procedimiento local requiere de más trabajo para alcanzar

1Este procedimiento adaptativo se conoce como global. Existe una alternativa local que procede subdividiendo los
intervalos de manera que cada uno de ellos tenga una estimación Ei del error acotada por |Ei| ≤ |(βi − αi)/(b− a)|

la tolerancia prefijada que el procedimiento global, pero, en compensasión, debería proveer, en general, resultados
que son aún más precisos.

2Mientras que una regla de Gauss de (2n + 1)-puntos es exacta para polinomios de hasta grado 4n + 1, una regla

2.

2

La implementación computacional del esquema de cuadratura adaptativo global con reglas de
Gauss-Kronrod es el núcleo del paquete de integración numérica conocido como quadpack, las
cuales a su vez forman parte del paquete de rutinas matemáticas slatec. Estas rutinas están
escritas en fortran 77 y entre ellas, la subrutina apropiada para nuestro problema planteado al
inicio de esta discusión, es la subrutina QAGS. La misma implementa un procedimiento adaptativo
y global de cuadratura en base una regla de Gauss-Kronrod de 10 y 21 puntos, respectivamente,
junto con un procedimiento de extrapolación. El uso de su implementación, en doble precisión,
DQAGS, es discutida a continuación.

Uso de la subrutina DQAGS del paquete QUADPACK

Para datos reales de doble precisión la interface de la subrutina DQAGS es como sigue:

SUBROUTINE DQAGS( F, A, B, EPSABS, EPSREL, RESULT, ABSERR, NEVAL,

IER, LIMIT, LENW, LAST, IWORK, WORK)

DOUBLE PRECISION
INTEGER
INTEGER
DOUBLE PRECISION

F, A, B, EPSABS, EPSREL, RESULT, ABSERR
NEVAL, IER, LIMIT, LENW, LAST
IWORK(LIMIT)
WORK(LENW)

donde,
F:

Es un dato de entrada que apunta a una función externa de doble precisión de una variable
de doble precisión, correspondiente al integrando f.
A:
Es un dato de entrada de doble precisión correspondiente al límite inferior a de la integral.
B:
Es un dato de entrada de doble precisión correspondiente al límite superior b de la integral.
EPSABS: Es un dato de entrada de doble precisión correspondiente a una tolerancia prefijada en

el error absoluto de la aproximación de la integral.

IER:

EPSREL: Es un dato de entrada de doble precisión correspondiente a una tolerancia prefijada en

el error relativo de la aproximación de la integral.

RESULT: Es un dato de salida de doble precisión que da la aproximación numérica Q de la integral.
ABSERR: Es un dato de salida de doble precisión que da una cota del error absoluto de Q, esto es,

|Q − I| ≤ ABSERR.

NEVAL: Es un dato de salida entero que da el número total de evaluaciones de f efectuadas por

el método.
Es un dato de salida entero utilizada como clave de éxito o error en el cómputo de la
aproximación numérica de la integral. Si IER=0, la subrutina ha calculado exitosamente
una aproximación numérica Q de la integral dentro de la tolerancia prefijada |Q − I| ≤
máx {EPSABS, EPSREL|I|} con una cota de error |Q−I| ≤ ABSERR. Valores de IERR no nulos
indican una terminación anormal del método y, en consecuencia, los valores devueltos
para la aproximación numérica y su estimación de error no son confiables. IERR = 1
significa que el númer