Clase No.

18:

Integración numérica: Regla del trapecio

Método de Romberg
MAT–251 Dr. Alonso Ramírez Manzanares
CIMAT A.C.
e-mail: alram@ cimat.mx
web: http://www.cimat.mx/salram/met_num/

Dr. Joaquín Peña Acevedo

CIMAT A.C.
e-mail: joaquin@ cimat.mx

Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 19.10.2015 1 / 19

Integración numérica

Dada una función f : [a, b] −→ R continua, queremos calcular la integral

definida de f en [a, b]:
Zb
f (x) dx
a

Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 19.10.2015 2 / 19

Integración numérica

Dada una función f : [a, b] −→ R continua, queremos calcular la integral

definida de f en [a, b]:
Zb
f (x) dx
a
Por el teorema fundamental del cálculo, si F es una primitiva de f , entonces
Zb
f (x) dx = F(b) − F(a).
a

Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 19.10.2015 2 / 19

Integración numérica

Dada una función f : [a, b] −→ R continua, queremos calcular la integral

definida de f en [a, b]:
Zb
f (x) dx
a
Por el teorema fundamental del cálculo, si F es una primitiva de f , entonces
Zb
f (x) dx = F(b) − F(a).
a

• No siempre se conoce una primitiva de f .

• De hecho, puede ser que la primitiva F no puede expresarse mediante
funciones elementales.
• La alternativa es realizar una estimación del valor de la integral.
• La integración numérica es conocida como cuadratura, que es el
nombre del método empírico en que se trazaban cuadrados bajo la
curva en cuestión para estimar el área.

Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 19.10.2015 2 / 19

Integración numérica

• En la práctica, sólo se usan algunos valores de la función f en puntos del

intervalo [a, b] para realizar la estimación de la integral.
• La evaluación de la expresión de F(x) puede ser menos estable
numéricamente que aplicar una regla de cuadratura de f (x).
• Aún cuando F(x) tenga una expresión algebraica estable conocida,
puede ser costosa evaluarla numéricamente.
Se puede decir que cambiamos el problema de integrar f por otro en el que
hay que calcular la integral de f + ∆f , que es más fácil de integral, eligiendo
∆f de modo que su integral no sea muy grande.

Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 19.10.2015 3 / 19

 Aproximación por sumas de Riemann (I)
Hacemos una partición P del intervalo

Pn : a = x0 < x1 < x2 < · · · < xn = b

Definimos
mi = min{f (x) : xi ≤ x ≤ xi+1 }
Mi = max{f (x) : xi ≤ x ≤ xi+1 }

n−1
X n−1
X
L(f , Pn ) = mi (xi+1 − xi ), U(f , Pn ) = Mi (xi+1 − xi )
i=0 i=0
2.5

2.5
2.0

2.0
1.5

1.5
c(−0.2, 2.5)

c(−0.2, 2.5)
1.0

1.0
0.5

0.5
0.0

0.0

a b a b

Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 19.10.2015 4 / 19

Aproximación por sumas de Riemann (II)
Se tiene que
• L(f , Pn ) ≤ U(f , Pn ).
• Si n aumenta, L(f , Pn ) debe aumentar, mientras que U(f , Pn ) debería
decrecer.
Tenemos que
Zb
L(f , Pn ) ≤ f (x) dx ≤ U(f , Pn ).
a
Para simplicar, podemos usar una partición uniforme de modo que
xi = a + b−a
n
i para i = 0, 1, ..., n.
Podemos aproximar el valor de la integral como
Zb
1
f (x) dx = [L(f , Pn ) + U(f , Pn )]
a 2
El error más grande que se comete con esta aproximación es
1
[U(f , Pn ) − L(f , Pn )]
2

Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 19.10.2015 5 / 19

Aproximación por sumas de Riemann (III)

El problema de usar este enfoque es que necesitamos calcular el máximo y

el mínimo de la función en cada subintervalo [xi , xi+1 ], lo cual no es práctico.

Ejemplo: Al calcular la integral de sin x de 0 a π/ 2 se obtienen los

siguientes resultados:
1
n L(f , Pn ) U(f , Pn ) 2
[L(f , Pn ) + U(f , Pn )]
100 0.99212546 1.00783342 0.99997944
1000 0.99921440 1.00078520 0.99999980
10000 0.99992146 1.00007854 1.00000000

Así, necesitamos otro tipo de aproximaciones de la integral:

• que no requieran hacer cálculos complicados sobre la función f ,
• que dependan de evaluar la función f en algunos puntos, y
• que ese número de puntos no sea demasiado grande.

Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 19.10.2015 6 / 19

Regla del trapecio (I)

La regla del trapecio aproxima la integral

Z xi+1
f (x) dx
xi

por el área del trapecio: 2.5

2.0
1.5
c(−0.2, 2.5)

1.0
0.5
0.0

xi xi+1
Z xi+1
0.0 0.5 1 1.0 1.5 2.0
Ii = f (x) dx ≈ [f (xi ) + f (xi+1 )](xi+1 − xi )
xi 2c(−0.2, 2.2)

Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 19.10.2015 7 / 19

Regla del trapecio (II)
Si tenemos la partición

b−a
a = x0 < x1 < x2 < · · · < xn = b, xi = a + i para i = 0, 1, ..., n.
n
b−a
Si definimos h = n
, entonces
Zb n−1
X n−1
X Z xi+1 h n−1
X
f (x) dx = Ii = f (x) dx ≈ [f (xi ) + f (xi+1 )]
a i=0 i=0 xi 2 i=0
Lo podemos reescribir como
 
Zb
f (a) + f (b) n−1
X
f (x) dx ≈ h  + f (xi )
a 2 i=1

Ejemplo n L(f , Pn ) U(f , Pn ) Trapecio

Aproximación 100 0.99212546 1.00783342 0.9999794382
de 1000 0.99921440 1.00078520 0.9999997944
R π/ 2
0
sin x dx 10000 0.99992146 1.00007854 0.9999999979

Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 19.10.2015 8 / 19

Regla del trapecio (III)

Para estimar el error que se comete con la regla del trapecio, sea pi (x) el
polinomio de grado a lo más 1 que interpola los puntos (xi , f (xi )) y
(xi+1 , f (xi+1 )). Entonces
Z xi+1
h h
Ti = pi (x) dx = [pi (xi+1 ) − pi (xi )] = [f (xi+1 ) − f (xi )]
xi 2 2
De lo visto en la parte de interpolación, el error entre pi (x) y f (x) para
x ∈ [xi , xi+1 ] es

1
f (x) − pi (x) = f 00 (ξx )(x − xi )(x − xi+1 ),
2
para algún ξx ∈ [xi , xi+1 ]
Z xi+1 Z xi+1
1
Ii − Ti = [f (x) − pi (x)] dx = f 00 (ξx )(x − xi )(x − xi+1 ) dx
xi 2 xi

Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 19.10.2015 9 / 19

Regla del trapecio (IV)

En cálculo, tenemos el siguiente resultado: Si α(x) es continua en [a, b] y

β(x) es integrable en [a, b] y no cambia de signo, entonces existe c ∈ [a, b]
tal que
Zb Zb
α(x)β(x) dx = α(c) β(x) dx
a a
En nuestro caso, (x − xi )(x − xi+1 ) es no positiva [xi , xi+1 ] en y podemos
suponer que f 00 es continua. Entonces
Z xi+1
h3
Ii − Ti = f 00 (ξi ) (x − xi )(x − xi+1 ) dx = − f 00 (ξi )
xi 12

para algún ξi ∈ [xi , xi+1 ].

Así, en todo el intervalo [a, b], el error es
 
n−1
X h3 n−1
X (b − a)h2 1 n−1
X
(Ii − Ti ) = − f 00 (ξi ) = −  f 00 (ξi ) .
i=0
12 i=0
12 n i=0

Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 19.10.2015 10 / 19

Regla del trapecio (V)
El promedio de las segundas derivadas es un valor que está dentro del
intervalo del valor mínimo y máximo de f 00 en [a, b]. Como supusimos que
f 00 es continua, debe existir un valor ξ ∈ [a, b] tal que

1 n−1
X
f 00 (ξ) = f 00 (ξi ).
n i=0

Así, el error cometido por la regla del trapecio es

(b − a)h2
− f 00 (ξ)
12

De este modo, si la discretización se hace más fina, el error de la

aproximación se reduce.

Para una discretización dada, podemos esperar que el error no sea muy
grande si el rango de valores de f 00 también es pequeño.

Ejemplo: Dar un valor para h, tal que el error cometido al estimar

R π/ 2
0
sin x dx con la regla de trapecio se menor que 10−6 .

Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 19.10.2015 11 / 19

Regla recursiva del trapecio (I)
Para una partición en el que el número de subintervalos es 2n , entonces
h = (b − a)/ 2n y podemos escribir la regla del trapecio de la siguiente forma:

b−a


 2 [f (a) + f (b)] n=0




R(n, 0) = n
2X −1
 h


2 [f (a) + f (b)] + h f (a + kh) n>0
k=1

Podemos dar una fórmula recursiva del método de trapecio:

n−1
1 2X
R(n, 0) = R(n − 1, 0) + h f (a + (2k − 1)h) n>0
2 k=1

Para ver esto, partimos de la identidad

1 1
 
R(n, 0) = R(n − 1, 0) + R(n, 0) − R(n − 1, 0)
2 2
y entonces calculamos la expresión que está entre corchetes.

Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 19.10.2015 12 / 19

Regla recursiva del trapecio (II)
h
Definimos C = 2
[f (a) + f (b)]. Entonces
n
2X −1
R(n, 0) = h f (a + ih) + C
i=1

2n−1
X−1
R(n − 1, 0) = 2h f (a + 2jh) + 2C
j=1

Entonces
n
1 2X −1 2n−1
X−1
R(n, 0) − R(n − 1, 0) = h f (a + ih) − h f (a + 2jh)
2 i=1 j=1
n−1
2X
= h f (a + (2k − 1)h)
k=1

La ventaja es que no tenemos que reevaluar el integrando en los puntos en

donde ya lo hemos evaluado.

Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 19.10.2015 13 / 19

Extrapolación de Richardson (I)

Consideramos una función escalar A(h) para la cual se cumple que existen
constantes a0 , a1 independientes de h, a1 6= 0, y dos números positivos p, p0 ,
con p0 > p, tales que
0
A(h) = a0 + a1 hp + O(hp ), h → 0.
Dado 0 < q < 1, entonces

0
A(q−1 h) = a0 + a1 q−p hp + O(hp )

Multiplicando la primer ecuación por q−p y restándole la segunda, se obtiene

0
q−p A(h) − A(q−1 h) = (q−p − 1)a0 + O(hp )

q−p A(h) − A(q−1 h) 0 q−p A(h) − A(h) + A(h) − A(q−1 h) 0

a0 = + O(hp ) = + O(hp )
q−p − 1 q−p − 1
A(h) − A(q−1 h) 0
a0 = A(h) + + O(hp ) h → 0.
q−p − 1

Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 19.10.2015 14 / 19

Extrapolación de Richardson (II)

Así, si tenemos dos aproximaciones de A(0), A(h) y A(q−1 h), que tienen el
mismo error O(hp , podemos obtener una mejor aproximación:

A(h) − A(q−1 h)
AER = A(h) +
q−p − 1
0
con un error menor de orden O(hp ), que es llamado extrapolación de
Richardson.

Esto ayuda a mejorar la aceleración en que converge una sucesión A(hn ) en

la que hn → 0, y A(0) es el valor que queremos obtener.

Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 19.10.2015 15 / 19

Método de Romberg (I)

Este método produce un arreglo triangular de la forma

R(0, 0)
R(1, 0) R(1, 1)
R(2, 0) R(2, 1) R(2, 2)
R(3, 0) R(3, 1) R(3, 2) R(3, 3)
.. .. .. .. ..
. . . . .
R(n, 0) R(n, 1) R(n, 2) R(n, 3) ··· R(n, n)

• La primera columna está dada por la aproximaciones de la integral dada

por la regla recursiva de trapecio.
• El resto de las columnas se obtienen por la fórmula de extrapolación de
Richardson

Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 19.10.2015 16 / 19

Método de Romberg (II)
Supongamos que tenemos

b−a
R(0, 0) = [f (a) + f (b)]
2
b−a a+b
   
R(1, 0) = f (a) + 2f + f (b)
4 2

Si I(f ) es el valor exacto de la integral, entonces

R(0, 0) = I(f ) + a1 h2 + O(h4 )

R(1, 0) = I(f ) + a1 (h/ 2)2 + O(h4 )

1
Mediante extrapolación de Richardson, con q = 2
y p = 2,

R(1, 0) − R(0, 0)
R(1, 1) = R(1, 0) + .
3

Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 19.10.2015 17 / 19

Método de Romberg (III)

Entonces I(f ) = R(1, 1) + O(h4 ). Supongamos que tenemos calculado R(2, 0).
Entonces podemos usar R(1, 0) y R(2, 0) para calcular R(2, 1). Así

R(1, 1) = I(f ) + a2 h4 + O(h6 )

R(2, 1) = I(f ) + a2 (h/ 2)4 + O(h6 )

Aplicando el método de extrapolación, podemos obtener la aproximación

R(2, 2).

De manera general, calculamos

1
R(n, m) = R(n, m − 1) + [R(n, m − 1) − R(n − 1, m − 1)]
4m −1
con n ≥ 1, m ≥ 1.

Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 19.10.2015 18 / 19

Método de Romberg (IV)

R π/ 2
Ejemplo: Estimar 0
sin x dx con el método de Romberg:

R(0, 0) = 0.78540
R(1, 0) = 0.94806 R(1, 1) = 1.002280
R(2, 0) = 0.98712 R(2, 1) = 1.000135 R(2, 2) = 0.999992
R(3, 0) = 0.99679 R(3, 1) = 1.000008 R(3, 2) = 0.999999 R(3, 3) = 1.000000008

Joaquín Peña (CIMAT) Métodos Numéricos (MAT–251) 19.10.2015 19 / 19