UNA INTRODUCCIÓN A LA CUADRATURA NUMÉRICA

illo
F. VADILLO

Resumen. Por razones históricas el cálculo aproximado de integrales se cono-

ce con el nombre de cuadratura mientras que el término integración se reserva
para la resolución de ecuaciones diferenciales. En esta breve introducción pri-
mero se presentan los métodos clásico: fórmulas del punto medio, trapecios

d
y regla de Simpson y se estudian los errores cometidos. Después se preset-
na las reglas de cuadratura compuestas, las fórmulas gaussianas y se naliza
estudiando los algoritmos de la cuadratura adaptativa.

1. Introducción Va
2.1. El método de interpolación
2.2. El método directo
2.3. El método de desarrollo de Taylor
Índice

2. Método habituales para obtener fórmulas de cuadratura

1
2
2
3
do
3
3. Error de cuadratura 4
3.1. El teorema del núcleo de Peano 5
4. Reglas de cuadratura compuestas 6
5. Cuadratura adaptativa 7
6. Cuadratura Gaussiana 8
an

6.1. Fórmulas de Randau y Lobatto 10

Referencias 10
rn

1. Introducción

En este capítulo se estudiarán métodos para el cálculo aproximado de integrales

denidas
Z b
(1.1) I(f ) = f (x) dx,
Fe

a
donde f es una función integrable en el intervalo [a, b]. Esta necesidad surge por
varias razones:
Sólo una pequeña parte de las funciones elementales tienen primitivas tam-
bién elementales. Considere por ejemplo la función exp(−x2 ).
Received by the editors 22 de febrero de 2021.
1
 2 F. VADILLO

Aunque exista la función primitiva, llega a ella puede resultar tedioso, pien-
se por ejemplo en las funciones racionales, o también pueden depender de
cambios de variables más o menos ingeniosos.
En muchas aplicaciones sólo se tiene una tabla de valores de f o una subru-
tina que permite evaluar f en puntos de su dominio.

illo
La mayoría de las aproximaciones numéricas llamadas fórmulas de cuadratura
son de la forma
n
(1.2)
X
In+1 (f ) = αj f (xj ),
j=0

donde los αj son los coecientes ó pesos y los puntos x0 < x1 < · · · < xn son los
nodos y el error de cuadratura es

d
(1.3) En+1 (f ) = I(f ) − In+1 (f ).
Para conocer la calidad de una fórmula de cuadratura se dene su grado de precisión
ó exactitud. Se dice que una fórmula de cuadratura tiene grado de exactitud M ≥ 0

2. Va
si calcula exactamente la integral para cada polinomio de grado menor o igual que
M , pero no es exacta para al menos un polinomio de grado M + 1.

Método habituales para obtener fórmulas de cuadratura

Suponiendo que los nodos x0 , · · · , xn están dados, por ejemplo si la función sólo
se conoce en una table. Para elegir los pesos se puede hacer de dos maneras
do
2.1. El método de interpolación. Se construye el polinomio de interpolación
pasando por lo nodos
n
X
Pn (x) = f (xj ) Lj (x),
j=0
y después
an

n Z b
(2.1)
X
In+1 (f ) = f (xj ) Lj (x) dx,
j=0 a

con los pesos αj = a Lj (x) dx para j = 0, · · · , n. Estas fórmulas de cuadratura se

Rb

las denomina formulas interpolatorias.

rn

Teorema 2.1. Dados los nodos x0 < x1 < · · · < xn , la fórmula interpolatoria tiene
un grado de exactitud mayor o igual que n.
Demostración. Si f es un polinomio de grado menor o igual que n, su polinomio
interpolador es también f y por tanto I(f ) = In+1 (f ). 
Fe

Ejemplo 2.2. Fórmula del rectángulo

IR (f ) = (b − a) f (a).
Ejemplo 2.3. Fórmula del punto medio
 
a+b
IP M (f ) = (b − a) f .
2
 UNA INTRODUCCIÓN A LA CUADRATURA NUMÉRICA 3

Ejemplo 2.4. Fórmula de los trapecios

b−a
IT (f ) = (f (a) + f (b)).
2

2.2. El método directo. Otra opción es determinar los pesos directamente es-
cribiendo el sistema de ecuaciones para los polinomios 1, x, x2 , · · · , xn que resulta

illo
α0 + α1 + · · · + αn = b − a,
b2 − a 2
α0 x0 + α1 x1 + · · · + αn xn = ,
2
.. ..
. = .
bn+1 − an+1
α0 xn0 + α1 xn1 + · · · + αn xnn = ,

d
n+1

un sistema con una única solución porque su matriz es de Valmermonde con abscisas
dos a dos distintas. De este modo se ha probado el siguiente teorema

Va
Teorema 2.5. Dados n + 1 nodos distintos dos a dos, existen unos únicos pesos
αj para que la regla (1.2) tenga un grado de exactitud mayo o igual a n
además,
Teorema 2.6. Dados n + 1 nodos distintos dos a dos, la regla (1.2) tiene un grado
de exactitud mayo o igual a n ssi es interpolatoria.
Es interesante destacar que en otras situaciones este método puede dar fórmulas
que no son interpolatorias, por ejemplo usando datos de la derivada de la función
do
en los nodos.
Ejemplo 2.7. Trate de encontrar una aproximación a f (x) dx con los datos
R1
−1
f (−1), f (0), f (1), f 00 (−1), f 00 (0) y f 00 (1).

2.3. El método de desarrollo de Taylor. En este tercer método se supone

an

que la función f es sucientemente diferenciable para desarrollar por Taylor hasta

donde sea necesario. Ilustraremos este método para obtener la regla de Simpsom:
 
a+b
IS (f ) = A f (a) + B f + C f (b).
2
Llamando c = a+b
2 yh= 2 ,
b−a
rn

b b
(x − c)2 00
Z Z 
0
I(f ) = f (x) dx = f (c) + (x − c)f (c) + f (c) + · · · dx
a a 2
h3 h5
= 2h f (c) + f 00 (c) + f (4) (c) + · · · .
3 60
Fe

Por otro lado

IS (f ) = A f (a) + B f (c) + C f (b)
h2
= (A + B + C) f (c) + h (−A + C) f 0 (c) + (A + C) f 00 (c) +
2
h3 h4
(−A + C) f 000 (c) + (A + C) f (4) (c) + . . . ,
6 24
 4 F. VADILLO

e identicando coecientes se tiene

2h = A + B + C,
0 = −A + C,
2
h = A + C,
3

illo
0 = −A + C,
2
h = A + C,
3
..
.
con la solución A = C = h
3 y B = 3 . Por tanto la regla de Simpson es:
4h
   
b−a a+b

d
IS (f ) = f (a) + 4 f + f (b) .
6 2

luego
3.

Va
Error de cuadratura

Suponiendo que se ha utilizado una fórmula interpolatoria para aproximar la

integral (1.1), el error cometido

En+1 (f ) = I(f ) − I(Pn ) =

Z

|En+1 (f )| ≤ (b − a) sup |f (x) − Pn (x)|,

a
b
(f (x) − Pn (x)) dx,
do
x∈[a,b]
de modo que si Pn es una buena aproximación a f , entonces In+1 (f ) es una buena
aproximación a I(f ).
Por otra parte, usando el error de interpolación f (x) − Pn (x) en la forma de
Newton se tiene
b b
f (n+1) (ξ(x))
Z Z
En+1 (f ) = f [x0 , · · · , xn , x] Π(x) dx =
an

Π(x) dx,
a a (n + 1)!
donde Π(x) = (x − x0 ) · · · (x − xn ). Entonces si Π(x) no cambia de signo en [a, b] y
f (n+1) es continua, se puede aplicar el teorema del valor medio del cálculo integral
para obtener
b
f (n+1) (α)
Z
En+1 (f ) = Π(x) dx,
rn

(n + 1)! a
con min{xj , a} < α < max{xj , b}.
Ejemplo 3.1. Para la fórmula del rectángulo
b
f 0 (α) (b − a)2
Z
(3.1) ER (f ) = (x − a) dx = f 0 (α) .
(1)! 2
Fe

Ejemplo 3.2. Para la regla de los trapecio

(b − a)3
(3.2) ET (f ) = −f 00 (α)
.
12
Cuando la función Π(x) cambia de signo en [a, b], este argumento no sirve y se
debe acudir a otros métodos si se desea expresar el error en términos de derivadas.
 UNA INTRODUCCIÓN A LA CUADRATURA NUMÉRICA 5

3.1. El teorema del núcleo de Peano. Peano descubrió un teorema que permi-
te expresar el error para una clase muy amplia de problemas numéricos. Se considera
el espacio C n+1 [a, b] donde se dene una forma lineal
Z b
L(f ) = [a0 (x)f (x) + a1 (x)f 0 (x) + · · · + an (x)f (n) (x)] dx

illo
a
n
X n
X n
X
+ bj,0 f (xj ) + bj,1 f 0 (xj ) + · · · + bj,n f (n) (xj ),
j=1 j=1 j=1

donde las funciones aj (x) son continuas a trozos y los nodos xj ∈ [a, b].
Si L(p) = 0 para todo polinomio en , entonces para cada f ∈
Q
Teorema 3.3. n
C n+1 [a, b]

d
Z b
L(f ) = f (n+1) (t) K(t) dt,
a
donde K(x) es el núcleo de Peano denido de la forma

Demostración. Ver [2, pag. 123].

1
n!

(x −
Va
Lx [ (x − t)n+ ],
K(t) =

y el subíndice denota que es aplicado como función de x de la forma

(x − t)n+
t)n+
=
=
(x − t)n ,
0,
x ≥ t,
x < t.

do
Corolario 3.4. Si además de las hipótesis del teorema anterior, el núcleo K(t) no
cambia de signo en [a, b] entonces para cada f ∈ C n+1 [a, b] existe un ξ ∈ [a, b] tal
que
f (n+1) (ξ)
L(f ) = L(xn+1 ).
(n + 1)!
an

Demostración. Por el teorema del valor medio

Z b
L(f ) = f (n+1) (ξ) K(t) dt,
a

que aplicado a xn+1 da

rn

Z b
L(xn+1 ) = (n + 1)! K(t) dt,
a

de donde se concluye de manera inmediata. 

Ejemplo 3.5. Se considera la regla del punto medio por lo que

Fe

Z b
a+b
EP M (f ) = f (x) dx − (b − a) f (c), con c = .
a 2
Su núcleo de Peanos es
Z b
K(t) = (x − t)+ dx − (b − a) (c − t)+ ,
a
 6 F. VADILLO

De modo que
1
Si t ≥ c ⇒ K(t) = (b − t)2 ,
2
1
Si t ≤ c ⇒ K(t) = (t − a)2 ,
2

illo
que al no cambiar el signo se puede aplicar el corolario y calculando
b
(a + b)2
Z
2 1
L(x ) = x2 dx − (b − a) = (b − a)3 ,
a 4 12
con lo que se llega a la expresión del error
f 00 (ξ)
(3.3) EP M (f ) = (b − a)3 .
24
Por un procedimiento parecido el error para la regla de Simpson es

d
Ejemplo 3.6.

f (4) (ξ)  b − a 5
(3.4) ES (f ) = − .
90 2

de forma que
4.

Z b
Va
Reglas de cuadratura compuestas

Para construir fórmulas de cuadratura compuesta se elige una partición del in-
tervalo de integración [a, b] de la forma
a = z0 < z1 < · · · < zn = b,

n−1 Z zj+1
do
X
f (x) dx = f (x) dx.
a j=0 zj

y después se aplica una cuadratura para cada integral del segundo miembro. Gene-
ralmente los sub intervalos son iguales de forma que zj = z0 + j h.
Aplicada por ejemplo la regla de los trapecios en cada sub intervalo resulta
h
an

[ f (z0 ) + 2 f (z1 ) + · · · + 2 f (zn−1 ) + f (zn ) ] .

2
Para conocer su error se usaría la expresión (??) con el resultado
n−1
h3 00
con zj ≤ ξj ≤ zj+1 .
X
ET (f ) = − f (ξj ),
j=0
12
rn

Lo que demostraría que si f 00 está acotada,

h3 h2
|ET (f )| ≤ n ||f 00 ||∞ = (b − a) ||f 00 ||∞ = O(h2 ).
12 12

De igual manera como

Fe

Z b n−1
X Z zj+1
f (x) dx = f (x) dx.
a j=1 zj−1

aplicando la regla de Simpson en cada subintervalo se tiene la regla de Simpson

compuesta
h
[ f (z0 ) + 4 f (z1 ) + 2 f (z2 ) + 4 f (z3 ) + · · · + 4 f (zn−1 ) + f (zn ) ] ,
3
 UNA INTRODUCCIÓN A LA CUADRATURA NUMÉRICA 7

con la acotación de error

b−a 4
|ES (f )| ≤ h ||f (4) ||∞ = O(h4 ).
90

illo
5. Cuadratura adaptativa

Considere por ejemplo, que necesite integrar la función f (x) = cos(x3 )200 en el
intervalo 0 ≤ x ≤ 3 representada en la Figura 1 con un error del orden de 10−6 . Si
utilizara la regla de Simpson y estimara el paso h en su expresión de error (3.4),
necesitaría del orden de 19200 sub-intervalos que evidentemente es una exageración

d
para calcular el área dibujada encerrada por la gráca de la función. Lo razonable
para este problema sería utilizar sub-intervalos pequeños en la zona de los picos y
sub-intervalos muchos mayores para las mesetas. Esta idea es la que desarrolla la

0.9

0.8

0.7

0.6

0.5

0.4
1

Va
técnicas de cuadratura adaptativas que se comentan.
do
0.3

0.2

0.1

0
0 0.5 1 1.5 2 2.5 3

Figura 1. Función cos(x3 )200 .

an

Suponga que quiere aproximar una integral denida I(f ) = f (x) dx utilizando
Rb
a
una fórmula de cuadratura Q(f ) tal que
(5.1) |Q(f ) − I(f )| < T ol,
rn

donde T ol representan una tolerancia dada.

Lo primero será realizar alguna estimación del error
E(f ; h) = K hq + O(hq+1 ),
suponiendo que q es el orden de la fórmula de cuadratura. Una técnica clásica es
evaluar la integral con dos pasos diferentes, por ejemplo h y h/2 y actuar de la
Fe

siguiente manera: si se utiliza una fórmula de los trapecios compuesta se tiene

h
R1 = [f (a) + 2f (a + h) + · · · + 2f (b − h) + f (b)] ,
2
h
R2 = [f (a) + 2f (a + h/2) + 2f (a + h) + · · · + 2f (b − h/2) + f (b)] ,
4
 8 F. VADILLO

con
I(f ) − R1 ≈ K h2 ,
 2
h
I(f ) − R2 ≈ K ,
2

illo
es decir
1
I(f ) − R2 ≈ ( I(f ) − R1 ).
4
por lo que
1
I(f ) − R1 = (I(f ) − R2 ) + (R2 − R1 ) ≈ ( I(f ) − R1 ) + (R2 − R1 ),
4
es decir

d
4
I(f ) − R1 ≈ ( R2 − R1 ),
3
1

Va
I(f ) − R2

I(f ) − S2

ver [4], [1] o también [6].

≈
que son las estimaciones a posteriori del error.

≈
3

15
1
15
( R2 − R1 ),

Para la regla de Simpson un análisis semejante daría:

I(f ) − S1
16
( S2 − S1 ),

( S2 − S1 )
do
Una vez que se tienen los estimadores a posteriori de los errores se actúa de la
siguiente manera: se comienza con una partición cualquiera del intervalo
a = x0 < x1 · · · < xr = b,
con hj = xj − xj−1 > 0 se computa Qj ≈ Ij = xjj−1 f (x) dx tal que
Rx
an

hj
|Ij − Qj | < T ol, j = 1, · · · , r,
b−a
esto se consiguePpartiendo el paso hasta cuando hiciera falta, y nalmente la apro-
ximación Qf = rj=1 Qj evidentemente vericará la condición (5.1).
El comando quad implementa en Matlab la regla de Simpson con paso adapta-
do, por ejemplo para función dibujada en la gura 1 el comando quad(@f,0,3,0.000001)
rn

devuelve la estimación 0,5316. Por otra parte quadgui de [5] muestra la evolución
del algoritmo y los sub-intervalos nales.

6. Cuadratura Gaussiana
Fe

Hasta ahora en las fórmulas de cuadratura (1.2) se suponía que los nodos xj es-
taban jados y lo que se podia decidir son los n + 1 pesos αj por lo que conseguían
fórmulas de exactitud n. Si también se pudieran elegir los n + 1 nodos x0 , · · · , xn ,
como se dispone de 2n + 2 parámetros, se buscarían fórmulas exactas para los po-
linomios 1, x, · · · , x2n+1 que son las denominas fórmulas de cuadratura Gaussianas
estudiadas en [7].
 UNA INTRODUCCIÓN A LA CUADRATURA NUMÉRICA 9

Evidentemente, los fórmulas Gaussianas son interpolatorias por lo que

Z b n
X Z b
f (x) dx = αj f (xj ) + f [x0 , · · · , xn , x] Πn+1 (x) dx,
a j=0 a

con Πn+1 (x) = (x − x0 ) · · · (x − xn ). Por tanto, los nodos deben vericar que

illo
Z b
(6.1) f [x0 , · · · , xn , x] Πn+1 (x) dx = 0,
a

para f (x) = 1, x, x2 , · · · x2n+1 .

Por otra parte, si p(x) ∈ Pm , es decir, es un polinomio de grado m, entonces
p(x) − p(x0 )
x 7→ p[x0 , x] = ∈ Pm−1 ,

d
x − x0
x 7→ p[x0 , x1 , x] ∈ Pm−2 ,
x 7→ p[x0 , x1 , x2 , x] ∈ Pm−3 ,

x
..
.

Va
7→ p[x0 , x1 , · · · , xn , x] ∈ Pm−(n+1) ,

Sustituyen entonces xn+1 , · · · , x2n+1 en (6.1) se concluye que Πn+1 (x) es un poli-
nomio ortogonal a todos los polinomios de grados menor o igual que n, es decir
Πn+1 (x) ⊥ Pn .
do
Además, (6.1) no puede ser cierta para polinomios de grado ≤ 2n + 2 pues entonces
Π(x) sería ortogonal a todo polinomio de grado ≤ n + 1 lo que es absurdo. En
conclusión:
Teorema 6.1. Una fórmula de cuadratura (1.2) usando n + 1 nodos puede tener
grado de precisión 2n + 1. Además, el grado máximo se alcanza únicamente cuando
an

usa fórmulas
Qn interpolatorias y los nodos x0 , x1 , · · · , xn son los ceros del polinomio
Πn+1 = j=0 (x − xj ) el n + 1-ésimo polinomio ortogonal en (a, b).

Además, si x0 , x1 , · · · , xn son los ceros de Πn+1 (x), entonces sus polinomios de

Lagrange son
Πn+1 (x)
rn

Lj (x) = ,
(x − xj ) Π0n+1 (xj )
y los pesos
Z b
1 Πn+1 (x)
(6.2) αj = dx,
Π0n+1 (xj ) (x − xj )
Fe

para j = 0, · · · , n.
Teorema 6.2. Los peso αj de la expresión (6.2) son positivos.
Demostración. Ver notas. 

En cuanto al error de cuadratura se demuestra el siguiente teorema:

 10 F. VADILLO

Teorema 6.3. Si la función f tiene 2n + 2 derivadas continuas en [a, b], el error

de cuadratura de la fórmula Gaussiana es
b
f (2n+2) (ξ)
Z
En+1 (f ) = Π2n+1 (x) dx.
(2n + 2)! a

illo
Demostración. Ver notas. 

6.1. Fórmulas de Randau y Lobatto. También se pueden considera fórmulas

de cuadratura con (n + 1) + (m + 1) del tipo
n m
(6.3)
X X
I(n+1)+(m+1) (f ) = αj f (xj ) + βj f (yj ),
j=0 j=0

d
con los nodos yj jos y los otros xj a elegir. En este caso, el orden que cabe esperar
es 2n + m + 1.
Los casos mas conocidos son
Fórmulas de Gauss: m = 0.

Va
Fórmula de Randau: m = 1, y1 = a.
Fórmula de Lobatto: m = 2, y1 = a, y2 = b.
Ver [3] para un detallado análisis.

Referencias
do
1. U.M. Ascher and C. Greif, A Firts Course in Numerical Methods, SIAM, 2011.
2. R. Bulirsch and J. Stoer, Introduction to Numerical Analysis, Springer, 1980.
3. P.J. Davis and P. Rabinowitz, Methods of Numerical Integration. Second Edition, Academic
Press, 1984.
4. M.H. Holmes, Introduction to Scientic Computing and Data Analysis, Springer, 2016.
5. C.B. Moler, Numerical Computing with Matlab, SIAM, 2004.
6. R. Sacco, F. Saleri, and A. Quarteroni, Numerical Mathematics, Springer, 2000.
an

7. J.M: Sanz Serna, La cuadratura gasussinan según Gauss, [Link]

pdf (2018).

Departamento de Matemáticas de la UPV/EHU

Email address : [Link]@[Link]
rn
Fe