CAPÍTULO

4 Diferenciación numérica e integración

Introducción
Una hoja de tejado corrugado se construye al presionar una hoja plana de aluminio dentro de
otra cuya sección transversal tiene la forma de una onda senoidal.

Se necesita una hoja corrugada de 4 pies de largo, la altura de cada onda es de 1 pulgada
desde la línea central y cada onda tiene un periodo de aproximadamente 2π. El problema de
encontrar la longitud de la hoja plana inicial consiste en encontrar la longitud de la curva de-
terminada por f (x) 5 sen x desde x 5 0 pulgadas hasta x 5 48 pulgadas. A partir del cálculo,
sabemos que esta longitud es
48 48
L= 1 + ( f (x))2 d x = 1 + (cos x)2 d x,
0 0

por lo que el problema se reduce a evaluar esta integral. A pesar de que la función senoidal
es una de las funciones matemáticas más comunes, el cálculo de su longitud implica una
integral elíptica de segunda clase, que no se puede evaluar de manera explícita. En este ca-
pítulo se desarrollan métodos para aproximar la solución a los problemas de este tipo. Este
problema particular se considera en el ejercicio 21 de la sección 4.4, en el ejercicio 15 de la
sección 4.5 y en el ejercicio 10 de la sección 4.7.
En la introducción del capítulo 3 mencionamos que una razón para usar polinomios
algebraicos para aproximar un conjunto arbitrario de datos es que, dada cualquier función
FRQWLQXDGHÀQLGDGHQWURGHXQLQWHUYDORFHUUDGRH[LVWHXQSROLQRPLRTXHHVWiDUELWUDULDPHQ-
te cerca de la función en cada punto del intervalo. Además, las derivadas y las integrales de
los polinomios se obtienen y evalúan con facilidad. No debería sorprender, entonces, que
muchos procedimientos para aproximar derivadas e integrales utilicen los polinomios
que aproximan la función.

127
128 CAPÍTULO 4 Diferenciación numérica e integración

4.1 Diferenciación numérica

La derivada de la función f en x0 es

f (x0 + h) − f (x0 )
f (x0 ) = lím .
h→0 h

Esta fórmula proporciona una forma obvia de generar una aproximación para f (x0 ); sim-
plemente calcule

f (x0 + h) − f (x0 )
h
para valores pequeños de h. Aunque esto puede ser obvio, no tiene mucho éxito debido a
nuestro antiguo némesis, el error de redondeo. Sin embargo, es un lugar para empezar.
Para aproximar f (x0 ), suponga primero que x0 ∈ (a, b), donde f ∈ C 2 [a, b], y
que x1 = x0 +h para alguna h = 0  TXH HV VXÀFLHQWHPHQWH SHTXHxD SDUD JDUDQWL]DU TXH
x1 ∈ [a, b]. Nosotros construimos el primer polinomio de Lagrange P0,1 (x) para f determi-
nado por x0 y x1, con su término de error:

(x − x0 )(x − x1 )
f (x) = P0,1 (x) + f (ξ(x))
2!
f (x0 )(x − x0 − h) f (x0 + h)(x − x0 ) (x − x0 )(x − x0 − h)
= + + f (ξ(x)),
−h h 2

para algunos ξ(x) entre x0 y x1. Derivando obtenemos

f (x0 + h) − f (x0 ) (x − x0 )(x − x0 − h)

f (x) = + Dx f (ξ(x))
h 2
f (x0 + h) − f (x0 ) 2(x − x0 ) − h
= + f (ξ(x))
h 2
(x − x0 )(x − x0 − h)
+ Dx ( f (ξ(x))).
2

Borrando los términos relacionados con ξ(x) obtenemos

f (x0 + h) − f (x0 )
f (x) ≈ .
h

8QDGLÀFXOWDGFRQHVWDIyUPXODHVTXHQRWHQHPRVLQIRUPDFLyQVREUH Dx f (ξ(x)), por lo

Isaac Newton usó y popularizó que el error de truncamiento no se puede calcular. Cuando x es x0VLQHPEDUJRHOFRHÀFLHQWH
las ecuaciones de diferencias de Dx f (ξ(x))HV\ODIyUPXODVHVLPSOLÀFDHQ
en el último cuarto del siglo
XVII, pero muchas de estas
f (x0 + h) − f (x0 ) h
técnicas fueron desarrolladas f (x0 ) = − f (ξ ). (4.1)
previamente por Thomas Harriot h 2
(1561–1621) y Henry Briggs
(1561–1630). Harriot realizó Para los valores pequeños de h, el cociente de diferencia [ f (x0 + h) − f (x0 )]/ h se
DYDQFHVVLJQLÀFDWLYRVHQWpFQLFDV puede utilizar para aproximar f (x0 ) con un error acotado por M|h|/2, donde M es una cota
de navegación y Briggs fue la
persona más responsable de
de | f (x)| para x entre x0 y x0 + h. A esta fórmula se le conoce como fórmula de diferencias
la aceptación de logaritmos como hacia adelante si h > YpDVHODÀJXUD \FRPRfórmula de diferencias hacia atrás
auxiliares para el cálculo. si h < 0.
 4.1 Diferenciación numérica 129

Figura 4.1
y
Pendiente f 9(x 0)
f (x0 1 h) 2 f (x 0)
Pendiente
h

x0 x0 1 h x

Ejemplo 1 Use la fórmula de diferencias hacia adelante para aproximar la derivada de f (x) = ln x
en x0 = 1.8 mediante h = 0.1, h = 0.05, y h = 0.01 y determine las cotas para los errores
de aproximación.

Solución La fórmula de diferencias hacia adelante

f (1.8 + h) − f (1.8)
h
con h = 0.1 nos da
ln 1.9 − ln 1.8 0.64185389 − 0.58778667
= = 0.5406722.
0.1 0.1
Puesto que f (x) = −1/x 2 y 1.8 < ξ < 1.9, una cota para este error de aproximación es
|h f (ξ )| |h| 0.1
= 2 < = 0.0154321.
2 2ξ 2(1.8)2

La aproximación y las cotas de error cuando h 5 0.05 y h 5 0.01 se encuentran de manera

similar y los resultados se muestran en la tabla 4.1.

Tabla 4.1 f (1.8 + h) − f (1.8) |h|

h f (1.8 + h)
h 2(1.8)2
0.1 0.64185389 0.5406722 0.0154321
0.05 0.61518564 0.5479795 0.0077160
0.01 0.59332685 0.5540180 0.0015432

Puesto que f (x) = 1/x , el valor exacto de f (1.8) es 0.555, y en este caso la cota del
error está bastante cerca del verdadero error de aproximación.

Para obtener las fórmulas generales de aproximación a la derivada, suponga que

{x0 , x1 , . . . , xn } son (n + 1) números distintos en algún intervalo I y que f ∈ C n+1 (I ). A
partir del teorema 3.3 en la página 83,
n
(x − x0 ) · · · (x − xn ) (n+1)
f (x) = f (xk )L k (x) + f (ξ(x)),
k=0
(n + 1)!
130 CAPÍTULO 4 Diferenciación numérica e integración

para algunos ξ(x) en I , donde L k (x) denota el k-ésimoFRHÀFLHQWHGHOSROLQRPLRGH/DJUDQ-

ge para f en x0 , x1 , . . . , xn. Al derivar esta ecuación obtenemos
n
(x − x0 ) · · · (x − xn )
f (x) = f (xk )L k (x) + Dx f (n+1) (ξ(x))
k=0
(n + 1!)
(x − x0 ) · · · (x − xn )
+ Dx [ f (n+1) (ξ(x))].
(n + 1)!

De nuevo tenemos un problema al calcular el error de truncamiento a menos que x sea

uno de los números xj. En este caso, el término que multiplica Dx [ f (n+1) (ξ(x))] es 0 y la
fórmula se vuelve
n n
f (n+1) (ξ(x j ))
f (x j ) = f (xk )L k (x j ) + (x j − xk ), (4.2)
k=0
(n + 1)! k=0
k= j

que recibe el nombre de fórmula de (n 1 1) puntos para aproximar f (x j ).

En general, el uso de más puntos de evaluación en la ecuación (4.2) produce mayor
precisión, a pesar de que el número de evaluaciones funcionales y el crecimiento del error de
redondeo disuade un poco esto. Las fórmulas más comunes son las relacionadas con tres y
cinco puntos de evaluación.
Primero derivamos algunas fórmulas útiles de tres puntos y consideramos aspectos de
sus errores. Puesto que
(x − x1 )(x − x2 ) 2x − x1 − x2
L 0 (x) = , tenemos L 0 (x) = .
(x0 − x1 )(x0 − x2 ) (x0 − x1 )(x0 − x2 )

De igual forma,
2x − x0 − x2 2x − x0 − x1
L 1 (x) = y L 2 (x) = .
(x1 − x0 )(x1 − x2 ) (x2 − x0 )(x2 − x1 )

Por lo tanto, a partir de la ecuación (4.2),

2x j − x1 − x2 2x j − x0 − x2
f (x j ) = f (x0 ) + f (x1 )
(x0 − x1 )(x0 − x2 ) (x1 − x0 )(x1 − x2 )
2
2x j − x0 − x1 1
+ f (x2 ) + f (3) (ξ j ) (x j − xk ), (4.3)
(x2 − x0 )(x2 − x1 ) 6 k=0
k= j

para cada j 5 0, 1, 2, donde la notación ξ j indica que este punto depende de xj.

Fórmulas de tres puntos

Las fórmulas a partir de la ecuación (4.3) se vuelven especialmente útiles si los nodos están
espaciados de manera uniforme, es decir, cuando
x1 = x0 + h y x2 = x0 + 2h, para algunas h = 0.
Supondremos nodos igualmente espaciados a lo largo del resto de esta sección.
Por medio de la ecuación (4.3) con x j = x0 , x1 = x0 + h, y x2 = x0 + 2h obtenemos
1 3 1 h 2 (3)
f (x0 ) = − f (x0 ) + 2 f (x1 ) − f (x2 ) + f (ξ0 ).
h 2 2 3

Al hacer lo mismo para x j = x1 obtenemos

1 1 1 h 2 (3)
f (x1 ) = − f (x0 ) + f (x2 ) − f (ξ1 )
h 2 2 6
 4.1 Diferenciación numérica 131

y, para x j = x2 ,

1 1 3 h 2 (3)
f (x2 ) = f (x0 ) − 2 f (x1 ) + f (x2 ) + f (ξ2 ).
h 2 2 3

Puesto que x1 = x0 + h y x2 = x0 + 2h, estas fórmulas también se pueden expresar

como

1 3 1 h 2 (3)
f (x0 ) = − f (x0 ) + 2 f (x0 + h) − f (x0 + 2h) + f (ξ0 ),
h 2 2 3
1 1 1 h 2 (3)
f (x0 + h) = − f (x0 ) + f (x0 + 2h) − f (ξ1 ), y
h 2 2 6
1 1 3 h 2 (3)
f (x0 + 2h) = f (x0 ) − 2 f (x0 + h) + f (x0 + 2h) + f (ξ2 ).
h 2 2 3

Por cuestiones de conveniencia, la sustitución de la variable x0 por x0 1 h se usa en

medio de la ecuación para cambiar esta fórmula por una aproximación para f 9(x0). Un cam-
bio similar, x0 para x0 1 2h, se utiliza en la última ecuación. Esto nos da tres fórmulas para
aproximar f 9(x0)

1 h 2 (3)
f (x0 ) = [−3 f (x0 ) + 4 f (x0 + h) − f (x0 + 2h)] + f (ξ0 ),
2h 3
1 h 2 (3)
f (x0 ) = [− f (x0 − h) + f (x0 + h)] − f (ξ1 ), y
2h 6
1 h 2 (3)
f (x0 ) = [ f (x0 − 2h) − 4 f (x0 − h) + 3 f (x0 )] + f (ξ2 ).
2h 3

Finalmente, observe que la última de estas ecuaciones se puede obtener a partir de la primera
simplemente al reemplazar h por 2h, de modo que en realidad sólo son dos fórmulas:

Fórmula del extremo de tres puntos

1 h 2 (3)
• f (x0 ) = [−3 f (x0 ) + 4 f (x0 + h) − f (x0 + 2h)] + f (ξ0 ) (4.4)
2h 3
donde j0 se encuentra entre x0 y x0 1 2h.

Fórmula del punto medio de tres puntos

1 h 2 (3)
• f (x0 ) = [ f (x0 + h) − f (x0 − h)] − f (ξ1 ), (4.5)
2h 6

donde j1 se encuentra entre x0 2 h y x0 1 h.

A pesar de que los errores en las dos ecuaciones (4.4) y (4.5) son O(h2), el error en la
ecuación (4.5) es aproximadamente la mitad del error en la ecuación (4.4). Esto porque
la ecuación (4.5) utiliza datos en ambos lados de x0 y la ecuación (4.4) los usa en un solo
lado. También observe que f se debe evaluar solamente en dos puntos en la ecuación (4.5),
PLHQWUDVTXHHQODHFXDFLyQ  VHQHFHVLWDQWUHVHYDOXDFLRQHV/DÀJXUDLOXVWUDODDSUR-
ximación producida a partir de la ecuación (4.5). La aproximación en la ecuación (4.4) es
útil cerca de los extremos de un intervalo porque la información sobre f fuera del intervalo
puede no estar disponible.
132 CAPÍTULO 4 Diferenciación numérica e integración

Figura 4.2
y

Pendiente f 9(x 0)

1
Pendiente [ f (x0 1 h) 2 f (x 0 2 h)]
2h

x0 2 h x0 x0 1 h x

Fórmulas de cinco puntos

Los métodos presentados en las ecuaciones (4.4) y (4.5) reciben el nombre de fórmulas de
tres puntos (aunque el tercer punto f (x0) no aparece en la ecuación (4.5)). De igual forma,
existen fórmulas de cinco puntos que implican la evaluación de la función en dos puntos
adicionales. El término de error para estas fórmulas es O(h4). Una fórmula de cinco pun-
tos común se usa para determinar las aproximaciones para la derivada en el punto medio.

Fórmula del punto medio de cinco puntos

1 h 4 (5)
• f (x0 ) = [ f (x0 − 2h) − 8 f (x0 − h) + 8 f (x0 + h) − f (x0 + 2h)] + f (ξ ),
12h 30
(4.6)
donde j se encuentra entre x0 2 2h y x0 1 2h.

La deducción de esta fórmula se considera en la sección 4.2. La otra fórmula de cinco puntos
se usa para aproximaciones en los extremos.

Fórmula del extremo de cinco puntos

1
• f (x0 ) = [−25 f (x0 ) + 48 f (x0 + h) − 36 f (x0 + 2h)
12h
h 4 (5)
+ 16 f (x0 + 3h) − 3 f (x0 + 4h)] + f (ξ ), (4.7)
5

donde j se encuentra entre x0 1 4h.

Las aproximaciones del extremo izquierdo se encuentran mediante esta fórmula con h > 0
y las aproximaciones del extremo derecho con h < 0. La fórmula del extremo de cinco puntos
es especialmente útil para la interpolación de spline cúbico condicionado de la sección 3.5.

Ejemplo 2 Los valores para f (x) = xe x se dan en la tabla 4.2. Utilice las fórmulas aplicables de tres y
cinco puntos para aproximar f (2.0).
 4.1 Diferenciación numérica 133

Solución Los datos en la tabla nos permiten encontrar cuatro diferentes aproximaciones de
Tabla 4.2
tres puntos. Podemos usar la fórmula del extremo (4.4) con h 5 0.1 o con h 5 20.1, y la
x f (x) fórmula del punto medio (4.5) con h 5 0.1 o con h 5 0.2.
1.8 10.889365 Mediante la fórmula del extremo (4.4) con h 5 0.1 obtenemos
1.9 12.703199
1
2.0 14.778112 [−3 f (2.0) + 4 f (2.1) − f (2.2] = 5[−3(14.778112) + 4(17.148957) − 19.855030)]
2.1 17.148957 0.2
2.2 19.855030 = 22.032310

y con h 5 20.1 nos da 22.054525.

Por medio de la fórmula del punto medio (4.5) con h 5 0.1 obtenemos

1
[ f (2.1) − f (1.9] = 5(17.148957 − 12.7703199) = 22.228790
0.2

y con h 5 0.2 nos da 22.414163.

/D~QLFDIyUPXODGHFLQFRSXQWRVSDUDODTXHODWDEODSURYHHGDWRVVXÀFLHQWHVHVODIyUPXOD
del punto medio (4.6) con h 5 0.1. Esto nos da

1 1
[ f (1.8) − 8 f (1.9) + 8 f (2.1) − f (2.2)] = [10.889365 − 8(12.703199)
1.2 1.2
+ 8(17.148957) − 19.855030]
= 22.166999.

Si no tenemos más información, aceptaríamos la aproximación del punto medio de cinco

puntos mediante h 5 0.1 como la más apropiada y esperamos que el valor verdadero se en-
cuentre entre esa aproximación y la aproximación del punto medio de tres puntos, es decir,
en el intervalo [22.166, 22.229].
En este caso, el valor verdadero es f (2.0) = (2 + 1)e2 = 22.167168, por lo que, en
realidad, los errores de aproximación son los siguientes:

Extremo de tres puntos con h = 0.1: 1.35 × 10−1 ;

Extremo de tres puntos con h = −0.1: 1.13 × 10−1 ;
Punto medio de tres puntos con h = 0.1: −6.16 × 10−2 ;
Punto medio de tres puntos con h = 0.2: −2.47 × 10−1 ;
Punto medio de cinco puntos con h = 0.1: 1.69 × 10−4 .

Los métodos también se pueden deducir para encontrar aproximaciones para derivadas
superiores de una función que sólo usa valores tabulados de la función en diferentes puntos.
La deducción es algebraicamente tediosa, sin embargo, sólo se presentará un procedimiento
representativo.
Represente una función f mediante la expansión de un tercer polinomio de Taylor sobre
un punto x0 y evalúe x0 + h y x0 2 h. Entonces,

1 1 1 (4)
f (x0 + h) = f (x0 ) + f (x0 )h + f (x0 )h 2 + f (x0 )h 3 + f (ξ1 )h 4
2 6 24
y

1 1 1 (4)
f (x0 − h) = f (x0 ) − f (x0 )h + f (x0 )h 2 − f (x0 )h 3 + f (ξ−1 )h 4 ,
2 6 24
donde x0 − h < ξ−1 < x0 < ξ1 < x0 + h.
134 CAPÍTULO 4 Diferenciación numérica e integración

Si adicionamos estas ecuaciones, los términos relacionados con f (x0 ) y − f (x0 ) se

cancelan, por lo que

1 (4)
f (x0 + h) + f (x0 − h) = 2 f (x0 ) + f (x0 )h 2 + [ f (ξ1 ) + f (4) (ξ−1 )]h 4 .
24

Al resolver esta ecuación para f (x0 ) obtenemos

1 h2
f (x0 ) = [ f (x0 − h) − 2 f (x0 ) + f (x0 + h)] − [ f (4) (ξ1 ) + f (4) (ξ−1 )]. (4.8)
h 2 24

Suponga que f (4) es continua en [x0 2 h, x0 1 h]. Puesto que 12 [ f (4) (ξ1 ) + f (4) (ξ−1 )] se
encuentra entre f (4) (ξ1 ) y f (4) (ξ−1 ), el teorema de valor intermedio implica que existe un
número j entre ξ1 y ξ−1 y, por lo tanto, en (x0 − h, x0 + h) con

1 (4)
f (4) (ξ ) = f (ξ1 ) + f (4) (ξ−1 ) .
2

(VWRQRVSHUPLWHUHHVFULELUODHFXDFLyQ  HQVXIRUPDÀQDO

Fórmula del punto medio de la segunda derivada

1 h 2 (4)
• f (x0 ) = [ f (x 0 − h) − 2 f (x 0 ) + f (x 0 + h)] − f (ξ), (4.9)
h2 12
Para algunos ξ , donde x0 − h < ξ < x0 + h.

Si f (4) es continua en [x0 2 h, x0 1 h], también está acotada, y la aproximación es O(h2).

Ejemplo 3 En el ejemplo 2 utilizamos los datos que se muestran en la tabla 4.3 para aproximar la pri-
mera derivada de f (x) 5 xex en x 5 2.0. Use la fórmula de la segunda derivada (4.9) para
Tabla 4.3 aproximar f 0(2.0).
x f (x)
Solución Los datos nos permiten determinar dos aproximaciones para f 0(2.0). Usando (4.9)
1.8 10.889365 con h 5 0.1 obtenemos
1.9 12.703199
2.0 14.778112 1
2.1 17.148957
[ f (1.9) − 2 f (2.0) + f (2.1)] = 100[12.703199 − 2(14.778112) + 17.148957]
0.01
2.2 19.855030
= 29.593200,

y mediante (4.9) con h 5 0.2 obtenemos

1
[ f (1.8) − 2 f (2.0) + f (2.2)] = 25[10.889365 − 2(14.778112) + 19.855030]
0.04
= 29.704275.

Puesto que f (x) = (x + 2)e x , el valor exacto es f (2.0) = 29.556224. Por lo tanto, los
errores reales son −3.70 × 10−2 y −1.48 × 10−1 , respectivamente.

Inestabilidad del error de redondeo

Es especialmente importante prestar atención al error de redondeo al aproximar derivadas.
Para ilustrar la situación, examinemos más de cerca la fórmula del punto medio de tres pun-
tos, la ecuación (4.5).

1 h 2 (3)
f (x0 ) = [ f (x0 + h) − f (x0 − h)] − f (ξ1 ),
2h 6
 4.1 Diferenciación numérica 135

Suponga que al evaluar f (x0 + h) y f (x0 − h), encontramos los errores de redondeo
e(x0 + h) y e(x0 − h). Entonces nuestros cálculos en realidad usan los valores f˜(x0 + h) y
f˜(x0 − h), que están relacionados con los valores verdaderos f(x0 1 h) y f(x0 2 h) mediante

f (x0 + h) = f˜(x0 + h) + e(x0 + h) y f (x0 − h) = f˜(x0 − h) + e(x0 − h).

El error total en la aproximación,

f˜(x0 + h) − f˜(x0 − h) e(x0 + h) − e(x0 − h) h 2 (3)

f (x0 ) − = − f (ξ1 ),
2h 2h 6

se debe tanto al error de redondeo, la primera parte, como al error de truncamiento. Si supo-
nemos que los errores de redondeo e(x0 ± h) están acotados por algún número ε > 0 y que
la tercera derivada de f está acotada por un número M > 0, entonces

f˜(x0 + h) − f˜(x0 − h) ε h2
f (x0 ) − ≤ + M.
2h h 6

Para reducir el error de truncamiento h 2 M/6, necesitamos reducir h. Pero conforme h se re-
duce, el error de redondeo ε/ h crece. En la práctica, entonces, casi nunca es ventajoso dejar
que h sea demasiado pequeña, porque en este caso, el error de redondeo dominará los cálculos.

Ilustración Considere utilizar los valores en la tabla 4.4. Para aproximar f (0.900), donde f (x) = sen x.
. El verdadero valor es cos 0.900 5 0.62161. La fórmula

f (0.900 + h) − f (0.900 − h)
f (0.900) ≈ ,
2h
con diferentes valores de h, proporciona las aproximaciones en la tabla 4.5.

Tabla 4.4 x sen x x sen x Tabla 4.5 Aproximación

h para f (0.900) Error
0.800 0.71736 0.901 0.78395
0.850 0.75128 0.902 0.78457 0.001 0.62500 0.00339
0.880 0.77074 0.905 0.78643 0.002 0.62250 0.00089
0.890 0.77707 0.910 0.78950 0.005 0.62200 0.00039
0.895 0.78021 0.920 0.79560 0.010 0.62150 −0.00011
0.898 0.78208 0.950 0.81342 0.020 0.62150 −0.00011
0.899 0.78270 1.000 0.84147 0.050 0.62140 −0.00021
0.100 0.62055 −0.00106

La mejor opción para h parece encontrarse entre 0.005 y 0.05. Podemos utilizar el cálcu-
ORSDUDYHULÀFDU FRQVXOWHHOHMHUFLFLR TXHVHSUHVHQWDXQPtQLPRSDUD

ε h2
e(h) = + M,
h 6
√
en h = 3
3ε/M, donde

M= máx | f (x)| = máx | cos x| = cos 0.8 ≈ 0.69671.

x∈[0.800,1.00] x∈[0.800,1.00]

Puesto que los valores de f están determinados para cinco lugares decimales, supondremos
que el error de redondeo está limitado por ε = 5 × 10−6. Por lo tanto, la mejor opción de h
es aproximadamente

3(0.000005)
h= ≈ 0.028,
3

0.69671
que es consistente con los resultados en la tabla 4.6.
136 CAPÍTULO 4 Diferenciación numérica e integración

En la práctica no podemos calcular una h óptima para utilizarla en la aproximación de

la derivada ya que no conocemos la tercera derivada de la función. Sin embargo, debemos
seguir siendo conscientes de que reducir el tamaño del paso no siempre mejora la aproxi-
mación.

Sólo hemos considerado los problemas de error de redondeo presentados por la fór-
PXODGHWUHVSXQWRVHFXDFLyQ  SHURVHSUHVHQWDQGLÀFXOWDGHVVLPLODUHVFRQWRGDVODV
fórmulas de diferenciación. La razón se puede rastrear a la necesidad de dividir entre una
potencia de h. Como encontramos en la sección 1.2 (consulte, especialmente, el ejemplo 3),
la división entre números pequeños tiende a exagerar el error de redondeo y, de ser posible,
debería evitarse esta operación. En caso de diferenciación numérica, no podemos evitar el
SUREOHPDSRUFRPSOHWRDSHVDUGHTXHORVPpWRGRVGHRUGHQVXSHULRUUHGXFHQODGLÀFXOWDG
Tome en cuenta que las Al igual que los métodos de aproximación, la diferenciación numérica es inestable ya
aproximaciones por método de
que los valores pequeños de h necesarios para reducir el error de truncamiento también
diferencias pueden ser inestables.
causan que el error de redondeo crezca. Esta es la primera clase de métodos inestables que
hemos encontrado y, de ser posible, estas técnicas deberían evitarse. Sin embargo, además
de usarse para propósitos computacionales, las fórmulas son necesarias para aproximar las
soluciones de ecuaciones ordinarias y diferenciales parciales.

La sección Conjunto de ejercicios 4.1 está disponible en línea. Encuentre la ruta de

acceso en las paginas preliminares

4.2 Extrapolación de Richardson

La extrapolación de Richardson se usa para generar resultados de alta precisión mientras usa
IyUPXODVGHEDMRRUGHQ$SHVDUGHTXHHOQRPEUHDGMXQWRDOPpWRGRVHUHÀHUHDXQDUWtFXOR
escrito por L. F. Richardson y J. A. Gaunt [RG] en 1927, la idea detrás de la técnica es mucho
más antigua. Un artículo interesante respecto a la historia y la aplicación de la extrapolación
se puede encontrar en [Joy].
La extrapolación se puede aplicar siempre que se sepa que una técnica de aproximación
tiene un término de error con un formato predecible, uno que depende de un parámetro, nor-
malmente el tamaño de paso h. Suponga que para cada número h = 0, tenemos una fórmula
N1 (h) que se aproxima a una constante M desconocida y que el error de truncamiento rela-
cionado con la aproximación tiene la forma

M − N1 (h) = K 1 h + K 2 h 2 + K 3 h 3 + · · · ,

para algún conjunto de constantes K 1 , K 2 , K 3 , . . . . (desconocidas).

Lewis Fry Richardson (1881- El error de truncamiento es O(h), a menos que haya una variación más grande entre las
1953) fue el primero en aplicar
sistemáticamente las matemáticas constantes K 1 , K 2 , K 3 , . . .,
a la predicción del tiempo
mientras trabajaba en Inglaterra M − N1 (0.1) ≈ 0.1K 1 , M − N1 (0.01) ≈ 0.01K 1 ,
SDUDOD2ÀFLQD0HWHRUROyJLFD
Como opositor concienzudo y, en general, M − N1 (h) ≈ K 1 h.
durante la Primera Guerra
Mundial, escribió ampliamente El objetivo de la extrapolación es encontrar una forma fácil de combinarlas en lugar de
sobre la futilidad económica de aproximaciones O(h) inapropiadas de manera adecuada para producir fórmulas con un error
la guerra, mediante sistemas de truncamiento de orden superior.
de ecuaciones diferenciales para Suponga, por ejemplo, que podemos combinar las fórmulas N1(h) para producir una
modelar interacciones racionales
entre países. La técnica de fórmula de aproximación O(h2), N2(h), para M con
extrapolación que lleva su
nombre fue el redescubrimiento M − N2 (h) = K̂ 2 h 2 + K̂ 3 h 3 + · · · ,
de una técnica con raíces que
son tan antiguas como Christiaan para algún, nuevamente desconocido, conjunto de constantes K̂ 2 , K̂ 3 , . . . . Entonces tenemos
Hugyens (1629–1695) y,
posiblemente, Arquímedes
(287–212 a.C.) M − N2 (0.1) ≈ 0.01 K̂ 2 , M − N2 (0.01) ≈ 0.0001 K̂ 2 ,
 4.2 Extrapolación de Richardson 137

y así sucesivamente. Si las constantes K 1 y K̂ 2 son aproximadamente de la misma magni-

tud, entonces las aproximaciones N2(h) serían mucho mejores que las aproximaciones N1(h)
correspondientes. La extrapolación continúa al combinar las aproximaciones N2(h) de mane-
ra que produce fórmulas con error de truncamiento O(h3) y así sucesivamente.
3DUD REVHUYDU HVSHFtÀFDPHQWH FyPR SRGHPRV JHQHUDU ODV IyUPXODV GH H[WUDSRODFLyQ
considere la fórmula O(h) para aproximar M
M = N1 (h) + K 1 h + K 2 h 2 + K 3 h 3 + · · · . (4.10)

Se asume que la fórmula se mantiene para todas las h positivas, por lo que reemplazamos el
parámetro h a la mitad de su valor. A continuación, tenemos una segunda fórmula de apro-
ximación O(h)
h h h2 h3
M = N1 + K1 + K2 + K3 + · · · . (4.11)
2 2 4 8

Restando dos veces la ecuación (4.11) de la ecuación (4.10) se elimina el término relaciona-
do con K1 y nos da

h h h2 h3
M = N1 + N1 − N2 (h) + K 2 − h2 + K3 − h3 + ··· .
2 2 2 4
(4.12)
'HÀQD

h h
N2 (h) = N1 + N1 − N1 (h) .
2 2

Entonces, la ecuación (4.12) es una fórmula de aproximación O(h2) para M:

K 2 2 3K 3 3
M = N2 (h) − h − h − ··· . (4.13)
2 4

Ejemplo 1 En el ejemplo 1 de la sección 4.1, utilizamos el método de diferencias hacia adelante con
h 5 0.1 y h 5 0.05 para encontrar aproximaciones para f (1.8) para f (x) 5 ln(x). Suponga
que esta fórmula tiene error de truncamiento O(h). Use la extrapolación en estos valores para
observar si esto resulta en una mejor aproximación.

Solución En el ejemplo 1 de la sección 4.1, encontramos que

con h = 0.1: f (1.8) ≈ 0.5406722, y con h = 0.05: f (1.8) ≈ 0.5479795.

Esto implica que

N1 (0.1) = 0.5406722 y N1 (0.05) = 0.5479795.

La extrapolación de estos resultados nos da la nueva aproximación

N2 (0.1) = N1 (0.05) + (N1 (0.05) − N1 (0.1)) = 0.5479795 + (0.5479795 − 0.5406722)

= 0.555287.

Se encontró que los resultados h 5 0.1 y h 5 0.05 son precisos dentro de 1.5 3 1022 y 7.7
3 1023, respectivamente. Puesto que f (1.8) = 1/1.8 = 0.5, el valor extrapolado es preciso
dentro de 2.7 3 1024.

La extrapolación se puede aplicar siempre que el error de truncamiento para una fórmula
tenga la forma
m−1
K j h α j + O(h αm )
j=1
138 CAPÍTULO 4 Diferenciación numérica e integración

para un conjunto de constantes Kj y cuando α1 < α2 < α3 < · · · < αm. Muchas fórmulas
utilizadas en la extrapolación tienen errores de truncamiento que sólo contienen potencias
pares de h, es decir, tienen la forma
M = N1 (h) + K 1 h 2 + K 2 h 4 + K 3 h 6 + · · · . (4.14)

La extrapolación es mucho más efectiva cuando todas las potencias de h están presentes de-
bido a que el proceso de promediar genera resultados con errores O(h 2 ), O(h 4 ), O(h 6 ), . . .,
con esencialmente, ningún incremento en el cálculo, sobre los resultados con errores
O(h), O(h 2 ), O(h 3 ), . . . .
Suponga que la aproximación tiene la forma de la ecuación (4.14). Al reemplazar h con
h / 2 obtenemos la fórmula de aproximación O(h2)

h h2 h4 h6
M = N1 + K1 + K2 + K3 + ··· .
2 4 16 64

Al restar cuatro veces esta ecuación de la ecuación (4.14), se elimina el término h2,

h h4 h6
3M = 4N1 − N1 (h) + K 2 − h4 + K3 − h6 + ··· .
2 4 16

Dividiendo esta ecuación entre 3, produce una fórmula O(h4)

1 h K2 h4 K3 h6
M= 4N1 − N1 (h) + − h4 + − h6 + ··· .
3 2 3 4 3 16

$OGHÀQLU
1 h h 1 h
N2 (h) = 4N1 − N1 (h) = N1 + N1 − N1 (h)
3 2 2 3 2

produce la fórmula de la aproximación con error de truncamiento O(h4):

h4 5h 6
M = N2 (h) − K 2 − K3 + ··· . (4.15)
4 16

Ahora reemplace h en la ecuación (4.15) con h/2 para producir una segunda fórmula O(h4)
h h4 5h 6
M = N2 − K2 − K3 − ··· .
2 64 1024

Al restar 16 veces esta ecuación de la ecuación (4.15) elimina el término h4 y da

h 15h 6
15M = 16N2 − N2 (h) + K 3 + ··· .
2 64

Al dividir esta ecuación entre 15 produce la nueva fórmula O(h6)

1 h h6
M= 16N2 − N2 (h) + K 3 + ··· .
15 2 64

Ahora tenemos la fórmula de aproximación O(h6)

1 h h 1 h
N3 (h) = 16N2 − N2 (h) = N2 + N2 − N2 (h) .
15 2 2 15 2
 4.2 Extrapolación de Richardson 139

Al continuar con este procedimiento obtenemos, para cada j = 2, 3, . . . , la aproximación

O(h2j)
h N j−1 (h/2) − N j−1 (h)
N j (h) = N j−1 + .
2 4 j−1 − 1

La tabla 4.6 muestra el orden en el que se generan las aproximaciones cuando

M = N1 (h) + K 1 h 2 + K 2 h 4 + K 3 h 6 + · · · . (4.16)

Se supone de manera conservadora que el verdadero resultado es preciso por lo menos

dentro del acuerdo de los dos resultados inferiores en la diagonal, en este caso, dentro de
|N3 (h) − N4 (h)|.

Tabla 4.6 O(h 2 ) O(h 4 ) O(h 6 ) O(h 8 )

1: N1 (h)
2: N1 ( h2 ) 3: N2 (h)
4: N1 ( h4 ) 5: N2 ( h2 ) 6: N3 (h)
7: N1 ( h8 ) 8: N2 ( h4 ) 9: N3 ( h2 ) 10: N4 (h)

Ejemplo 2 El teorema de Taylor se puede utilizar para mostrar que la fórmula de diferencias centradas
en la ecuación (4.5) para aproximar f (x0 ) se puede expresar con una fórmula de error:

1 h2 h 4 (5)
f (x0 ) = [ f (x0 + h) − f (x0 − h)] − f (x0 ) − f (x0 ) − · · · .
2h 6 120

Encuentre las aproximaciones de orden O(h2), O(h4), y O(h6) para f 9(2.0) cuando f 5 xex y
h 5 0.2.

Solución Probablemente, las constantes K 1 = − f (x0 )/6, K 2 = − f (5) (x0 )/120, · · · , se-
rán desconocidas, sin embargo, esto no es importante. Sólo necesitamos saber que estas
FRQVWDQWHVH[LVWHQFRQHOÀQGHDSOLFDUODH[WUDSRODFLyQ
Tenemos la aproximación O(h2)

h2 h 4 (5)
f (x0 ) = N1 (h) − f (x0 ) − f (x0 ) − · · · , (4.17)
6 120

donde

1
N1 (h) = [ f (x0 + h) − f (x0 − h)].
2h
Esto nos da las primeras aproximaciones O(h2)

1
N1 (0.2) = [ f (2.2) − f (1.8)] = 2.5(19.855030 − 10.889365) = 22.414160
0.4

1
N1 (0.1) = [ f (2.1) − f (1.9)] = 5(17.148957 − 12.703199) = 22.228786.
0.2
Al combinarlas para producir la primera aproximación O(h4) obtenemos

1 1
N2 (0.2) = N1 (0.1) + (N1 (0.1) − N1 (0.2)) = 22.228786 + (22.228786 − 22.414160)
3 3
= 22.166995.
 4.2 Extrapolación de Richardson 141

1 1
f (x0 − h) = f (x0 ) − f (x0 )h + f (x0 )h 2 − f (x0 )h 3
2 6
1 (4) 1
+ f (x0 )h 4 − f (5) (ξ2 )h 5 , (4.19)
24 120

donde x0 − h < ξ2 < x0 < ξ1 < x0 + h.

Al restar la ecuación (4.19) de la ecuación (4.18) obtenemos una nueva aproximación
para f (x)):

h3 h 5 (5)
f (x0 + h) − f (x0 − h) = 2h f (x0 ) + f (x0 ) + [ f (ξ1 ) + f (5) (ξ2 )], (4.20)
3 120
lo cual implica que
1 h2 h 4 (5)
f (x0 ) = [ f (x0 + h) − f (x0 − h)] − f (x0 ) − [ f (ξ1 ) + f (5) (ξ2 )].
2h 6 240
Si f (5) es continua en [ x0 − h, x0 + h], el teorema del valor intermedio 1.11 implica que
existe un número ξ̃ en (x0 − h, x0 + h) con
1 (5)
f (5) (ξ̃ ) = f (ξ1 ) + f (5) (ξ2 ) .
2
Como consecuencia, tenemos la aproximación O(h2)
1 h2 h 4 (5)
f (x0 ) = [ f (x0 + h) − f (x0 − h)] − f (x0 ) − f (ξ̃ ). (4.21)
2h 6 120
A pesar de que la aproximación en la ecuación (4.21) es igual a la que se obtuvo en la
fórmula de tres puntos en la ecuación (4.5), ahora, el punto desconocido de evaluación se
presenta en f (5) en lugar de en f -. La extrapolación aprovecha esto al reemplazar primero h
en la ecuación (4.21) con 2h para producir la nueva fórmula
1 4h 2 16h 4 (5)
f (x0 ) = [ f (x0 + 2h) − f (x0 − 2h)] − f (x0 ) − f (ξ̂ ), (4.22)
4h 6 120
donde ξ̂ se encuentra entre x0 − 2h y x0 + 2h.
Al multiplicar la ecuación (4.21) por 4 y restar la ecuación (4.22) produce
2 1
3 f (x0 ) = [ f (x0 + h) − f (x0 − h)] − [ f (x0 + 2h) − f (x0 − 2h)]
h 4h
h 4 (5) 2h 4 (5)
− f (ξ̃ ) + f (ξ̂ ).
30 15
Incluso si f (5) es continua en [x0 − 2h, x0 + 2h], el teorema de valor intermedio 1.11 no
se puede aplicar como lo hicimos para derivar la ecuación (4.21) porque tenemos la diferen-
cia de términos relacionados con f (5). Sin embargo, es posible usar un método alterno para
mostrar que f (5) (ξ̃ ) y f (5) (ξ̂ ) se puede seguir reemplazando con un valor común f (5) (ξ ).
Al suponer esto y dividir entre 3 produce la fórmula del punto medio de cinco puntos la
ecuación (4.6) que observamos en la sección 4.1:

1 h 4 (5)
f (x0 ) = [ f (x0 − 2h) − 8 f (x0 − h) + 8 f (x0 + h) − f (x0 + 2h)] + f (ξ ).
12h 30

Otras fórmulas para las primeras derivadas y las derivadas superiores se pueden deducir
de manera similar. Consulte, por ejemplo, el ejercicio 8.
A lo largo del texto se utiliza la técnica de extrapolación. Las aplicaciones más promi-
nentes se presentan al aproximar las integrales en la sección 4.5 y al determinar soluciones
aproximadas para las ecuaciones diferenciales en la sección 5.8.
142 CAPÍTULO 4 Diferenciación numérica e integración

La sección Conjunto de ejercicios 4.2 está disponible en línea. Encuentre la ruta de

acceso en las páginas preliminares.

4.3 Elementos de integración numérica

$PHQXGRVXUJHODQHFHVLGDGGHHYDOXDUODLQWHJUDOGHÀQLGDGHXQDIXQFLyQTXHQRWLHQHXQD
antiderivada o cuya antiderivada no es fácil de obtener. El método básico asociado con la
b
aproximación de a f (x) d x recibe el nombre de cuadratura numérica. Éste utiliza una
n b
suma i=0 ai f (xi ) para aproximar a f (x) d x.
Los métodos de cuadratura en esta sección se basan en los polinomios de interpolación
que se han explicado en el capítulo 3. La idea básica es seleccionar un conjunto de nodos
distintos {x0 , . . . , xn } del intervalo [a, b]. Entonces integramos el polinomio interpolante de
Lagrange
n
Pn (x) = f (xi )L i (x)
i=0

y su término de error de truncamiento sobre [a, b] para obtener

b b n b n
f (n+1) (ξ(x))
f (x) d x = f (xi )L i (x) d x + (x − xi ) dx
a a i=0 a i=0
(n + 1)!
n b n
1
= ai f (xi ) + (x − xi ) f (n+1) (ξ(x)) d x,
i=0
(n + 1)! a i=0

donde j(x) se encuentra en [a, b] para cada x y

b
ai = L i (x) d x, para cada i = 0, 1, . . . , n.
a

La fórmula de cuadratura es, por lo tanto,

b n
f (x) d x ≈ ai f (xi ),
a i=0

con un error dado por

b n
1
E( f ) = (x − xi ) f (n+1) (ξ(x)) d x.
(n + 1)! a i=0

Antes de analizar la situación general de las fórmulas de cuadratura, consideremos las

fórmulas producidas mediante el uso del primer y del segundo polinomios de Lagrange con
nodos igualmente espaciados. Esto da la regla trapezoidal y la regla de Simpson, las cuales
se presentan generalmente en cursos de cálculo.

La regla trapezoidal
b
Para derivar la regla trapezoidal (o regla del trapecio) para aproximar a f (x) d x, sean
x0 = a, x1 = b, h = b − a y utilice el polinomio de Lagrange

(x − x1 ) (x − x0 )
P1 (x) = f (x0 ) + f (x1 ).
(x0 − x1 ) (x1 − x0 )
 4.3 Elementos de integración numérica 143

Entonces
b x1
(x − x1 ) (x − x0 )
f (x) d x = f (x0 ) + f (x1 ) d x
a x0 (x0 − x1 ) (x1 − x0 )
(4.23)
x1
1
+ f (ξ(x))(x − x0 )(x − x1 ) d x.
2 x0

El producto (x − x0 )(x − x1 ) no cambia de signo en [x0, x1], por lo que el teorema del valor
promedio ponderado para integrales 1.13 se puede aplicar al término de error para obtener,
para algunos j en (x0, x1),
x1 x1
f (ξ(x))(x − x0 )(x − x1 ) d x = f (ξ ) (x − x0 )(x − x1 ) d x
x0 x0
Cuando usamos el término x1
trapezoidal, nos referimos a x3 (x1 + x0 ) 2
= f (ξ ) − x + x0 x1 x
XQDÀJXUDGHFXDWURODGRVFRQ 3 2 x0
al menos dos lados paralelos. El
WpUPLQRHXURSHRSDUDHVWDÀJXUD h3
es trapezium. Para confundir =− f (ξ ).
más el tema, la palabra europea 6
trapezoidalVHUHÀHUHDXQDÀJXUD
de cuatro lados sin ningún lado Por consiguiente, la ecuación (4.23) implica que
igual y la palabra estadounidense
b x1
SDUDHVWHWLSRGHÀJXUDHV (x − x1 )2 (x − x0 )2 h3
trapezium. f (x) d x = f (x0 ) + f (x1 ) − f (ξ )
a 2(x0 − x1 ) 2(x1 − x0 ) x0 12
(x1 − x0 ) h 3
= [ f (x0 ) + f (x1 )] − f (ξ ).
2 12

Por medio de la notación h = x1 − x0 obtenemos la siguiente regla:

Regla trapezoidal:
b
h h3
f (x) d x =
[ f (x0 ) + f (x1 )] − f (ξ ).
a 2 12
Esto recibe el nombre de regla trapezoidal porque cuando f es una función con valores
b
positivos, a f (x) d x se aproxima mediante el área de un trapecio, como se muestra en la
ÀJXUD

Figura 4.3
y
y 5 f (x)

y 5 P1(x)

a 5 x0 x1 5 b x

El término de error para la regla trapezoidal implica f 0, por lo que la regla da el resultado
exacto cuando se aplica a cualquier función cuya segunda derivada es idénticamente cero, es
decir, cualquier polinomio de grado uno o menos.
144 CAPÍTULO 4 Diferenciación numérica e integración

Regla de Simpson
La regla de Simpson resulta de la integración sobre [a, b] del segundo polinomio de
Lagrange con nodos igualmente espaciados x0 = a, x2 = b, y x1 = a + h, en don-
de h = (b − a)/2. (véase ODÀJXUD 

Figura 4.4
y
y 5 f (x)

y 5 P2(x)

a 5 x0 x1 x2 5 b x

Por lo tanto,
b x2
(x − x1 )(x − x2 ) (x − x0 )(x − x2 )
f (x) d x = f (x0 ) + f (x1 )
a x0 (x0 − x1 )(x0 − x2 ) (x1 − x0 )(x1 − x2 )
(x − x0 )(x − x1 )
+ f (x2 ) d x
(x2 − x0 )(x2 − x1 )
x2
(x − x0 )(x − x1 )(x − x2 ) (3)
+ f (ξ(x)) d x.
x0 6

Al deducir la regla de Simpson de esta forma, sin embargo, da un solo término de error O(h4)
relacionado con f (3). Al aproximar el problema de otra forma, se puede derivar otro término
de orden superior relacionado con f (4).
Para ilustrar este método alternativo, suponga que f se expande en el tercer polinomio de
Taylor alrededor de x1. Entonces, para cada x en [x0, x2], existe un número j(x) en (x0, x2) con

f (x1 ) f (x1 )
f (x) = f (x1 ) + f (x1 )(x − x1 ) + (x − x1 )2 + (x − x1 )3
2 6
f (4) (ξ(x))
+ (x − x1 )4
24
y
x2
f (x1 ) f (x1 )
f (x) d x = f (x1 )(x − x1 ) + (x − x1 )2 + (x − x1 )3
x0 2 6
x2 x2
f (x1 ) 1
+ (x − x1 )4 + f (4) (ξ(x))(x − x1 )4 d x. (4.24)
24 x0 24 x0

Puesto que (x 2 x1)4 nunca es negativo en [x0, x2], el teorema de valor promedio ponderado
para las integrales 1.13 implica que
x2 x2 x2
1 f (4) (ξ1 ) f (4) (ξ1 )
f (4) (ξ(x))(x − x1 )4 d x = (x − x1 )4 d x = (x − x1 )5 ,
24 x0 24 x0 120 x0
 4.3 Elementos de integración numérica 145

para algún número j1 en (x0, x2).

Sin embargo, h 5 x2 2 x1 5 x1 2 x0, por lo que

(x2 − x1 )2 − (x0 − x1 )2 = (x2 − x1 )4 − (x0 − x1 )4 = 0,

mientras

(x2 − x1 )3 − (x0 − x1 )3 = 2h 3 y (x2 − x1 )5 − (x0 − x1 )5 = 2h 5 .

Por consiguiente, la ecuación (4.24) se puede reescribir como

x2
h3 f (4) (ξ1 ) 5
f (x) d x = 2h f (x1 ) + f (x1 ) + h .
x0 3 60

Ahora, si reemplazamos f (x1 ) por medio de la aproximación determinada en la ecua-

ción (4.9) de la sección 4.1, tenemos
x2
h3 1 h 2 (4) f (4) (ξ1 ) 5
f (x) d x = 2h f (x1 ) + [ f (x0 ) − 2 f (x1 ) + f (x2 )] − f (ξ2 ) + h
x0 3 h 2 12 60
h h 5 1 (4) 1
= [ f (x0 ) + 4 f (x1 ) + f (x2 )] − f (ξ2 ) − f (4) (ξ1 ) .
3 12 3 5

Con métodos alternos se puede mostrar (consulte el ejercicio 26) que los valores ξ1 y ξ2 en
Thomas Simpson (1710–1761) esta expresión se pueden reemplazar mediante un valor común ξ en (x0 , x2 ). Esto da la regla
fue un matemático autodidacta de Simpson.
que en sus primeros años se
ganaba la vida como tejedor. Su Regla de Simpson:
principal interés fue la teoría
x2
de la probabilidad, aunque en h h 5 (4)
1750 publicó un libro de cálculo f (x) d x = [ f (x0 ) + 4 f (x1 ) + f (x2 )] − f (ξ ).
x0 3 90
de dos volúmenes titulado
The Doctrine and Application
of Fluxions (La doctrina y El término de error en la regla de Simpson implica la cuarta derivada de f, por lo que da
DSOLFDFLyQGHÁX[LRQHV . resultados exactos cuando se aplica a cualquier polinomio de grado tres o menos.
2
Ejemplo 1 Compare las aproximaciones de la regla trapezoidal y de la regla de Simpson para 0 f (x) d x
cuando f (x) es

a) x2
√ b) x4 c) (x + 1)−1
d) 1 + x2 e) sen x f) ex

Solución En [0, 2], las reglas trapezoidal y de Simpson tiene las formas
2
Trapezoidal: f (x) d x ≈ f (0) + f (2) y
0
2
1
De Simpson: f (x) d x ≈ [ f (0) + 4 f (1) + f (2)].
0 3

Cuando f (x) = x 2, obtenemos

2
Trapezoidal: f (x) d x ≈ 02 + 22 = 4 y
0
2
1 2 8
De Simpson: f (x) d x ≈ [(0 ) + 4 · 12 + 22 ] = .
0 3 3
La aproximación a partir de la regla de Simpson es exacta porque su error de truncamiento
implica f (4), lo cual es idénticamente 0 cuando f (x) = x 2 .
Los resultados a los tres lugares para las funciones se resumen en la tabla 4.7. Observe
TXHHQFDGDLQVWDQFLDODUHJODGH6LPSVRQHVVLJQLÀFDWLYDPHQWHVXSHULRU
146 CAPÍTULO 4 Diferenciación numérica e integración

Tabla 4.7 (a) (b) (c) (d) (e) (f)

√
f (x) x2 x4 (x + 1)−1 1 + x2 sen x ex
Valor exacto 2.667 6.400 1.099 2.958 1.416 6.389
Trapezoidal 4.000 16.000 1.333 3.326 0.909 8.389
De Simpson 2.667 6.667 1.111 2.964 1.425 6.421

Precisión de medición
Las derivadas estándar de las fórmulas de error de cuadratura están basadas al determinar la
clase de polinomios para los que estas fórmulas producen resultados exactos. La siguiente
GHÀQLFLyQVHXWLOL]DSDUDIDFLOLWDUHODQiOLVLVGHHVWDGHULYDGD
Definición 4.1 El grado de precisión, o precisión, de una fórmula de cuadratura es el mayor entero positivo n,
de tal forma que la fórmula es exacta para xk, para cada k = 0, 1, . . . , n.
/DGHÀQLFLyQLPSOLFDTXHODVUHJODVWUDSH]RLGDO\GH6LPSVRQWLHQHQJUDGRVGHSUH-
La precisión mejorada de la
cisión uno y tres, respectivamente.
regla de Simpson sobre la regla
trapezoidal se explica de manera La integración y la sumatoria son operaciones lineales; es decir,
intuitiva con el hecho de que la b b b
regla de Simpson incluye una
evaluación de punto medio que
(α f (x) + βg(x)) d x = α f (x) d x + β g(x) d x
a a a
proporciona mejor equilibrio
para la aproximación.
y
n n n
(α f (xi ) + βg(xi )) = α f (xi ) + β g(xi ),
i=0 i=0 i=0

para cada par de funciones integrables f y g y cada par de constantes reales a y b. Esto
implica (consulte el ejercicio 25) que

• el grado de precisión de una fórmula de cuadratura es n si y sólo si el error es cero para

La terminología abierta y cerrada
todos los polinomios de grado k = 0, 1, . . . , n, pero no es cero para algunos polinomios
para los métodos implica que el
método abierto sólo utiliza como de grado n 1 1.
nodos los puntos en el intervalo
abierto (a, b) para aproximar Las reglas trapezoidal y de Simpson son ejemplos de una clase de métodos conocidos como
b
a f (x) d x. Los métodos fórmulas de Newton-Cotes. Existen dos tipos de fórmulas de Newton-Cotes: abiertas y ce-
cerrados incluyen los puntos a rradas.
y b del intervalo cerrado [a, b]
como nodos.
Fórmulas de Newton-Cotes cerradas
La fórmula cerrada de (n 1 1) puntos de Newton-Cotes utiliza nodos xi = x0 +i h, para
i = 0, 1, . . . , n, donde x0 = a, xn = b y h = (b − a)/n.  9pDVH OD ÀJXUD   5HFLEH HO
nombre de cerrada porque los extremos del intervalo cerrado [a, b] se incluyen como nodos.

Figura 4.5 y

y = Pn(x)
y = f (x)

a 5 x0 x1 x2 xn21 xn 5 b x
 4.3 Elementos de integración numérica 147

La fórmula asume la forma

b n
f (x) d x ≈ ai f (xi ),
a i=0

donde
xn xn n
(x − x j )
ai = L i (x) d x = d x.
x0 x0 j=0
(xi − x j )
j =i

El siguiente teorema detalla el análisis de error relacionado con las fórmulas de

Newton-Cotes cerradas. Para una demostración de este teorema, consulte [IK], p. 313.
n
Teorema 4.2 Suponga que i=0 ai f (xi ) denota la fórmula cerrada de (n 1 1) puntos de Newton-Cotes
con x0 = a, xn = b, y h = b − a)/n. Existe ξ ∈ (a, b) para el que
Roger Cotes (1682–1716)
ascendió desde un origen b n n
h n+3 f (n+2) (ξ )
humilde hasta convertirse, en f (x) d x = ai f (xi ) + t 2 (t − 1) · · · (t − n) dt,
1704, en el primer profesor a i=0
(n + 2)! 0
plumiano en la Universidad
de Cambridge. Realizó
numerosos avances en las áreas si n es par y f ∈ C n+2 [a, b], y
de matemáticas, incluyendo
b n n
los métodos numéricos para h n+2 f (n+1) (ξ )
la interpolación e integración. f (x) d x = ai f (xi ) + t (t − 1) · · · (t − n) dt,
Newton es famoso por decir a i=0
(n + 1)! 0
respecto a Cotes: “Si hubiera
vivido, habríamos aprendido
algo”. si n es impar y f ∈ C n+1 [a, b].

Observe que cuando n es un entero par, el grado de precisión es n 1 1, a pesar de que el

polinomio de interpolación es de grado a lo sumo n. Cuando n es impar, el grado de precisión
sólo es n.
Se listan algunas de las fórmulas comunes de Newton-Cotes cerradas con sus términos
de error. Observe que, en cada caso, el valor desconocido j se encuentra en (a, b).

n 5 1: Regla trapezoidal

x1
h h3
f (x) d x = [ f (x0 ) + f (x1 )] − f (ξ ), donde x0 < ξ < x1 . (4.25)
x0 2 12

n 5 2: Regla de Simpson

x2
h h 5 (4)
f (x) d x = [ f (x0 ) + 4 f (x1 ) + f (x2 )] − f (ξ ), donde x0 < ξ < x2 .
x0 3 90
(4.26)

n 5 3: Regla de tres octavos de Simpson

x3
3h 3h 5 (4)
f (x) d x = [ f (x0 ) + 3 f (x1 ) + 3 f (x2 ) + f (x3 )] − f (ξ ), (4.27)
x0 8 80
donde x0 < ξ < x3 .
148 CAPÍTULO 4 Diferenciación numérica e integración

n 5 4:

x4
2h 8h 7 (6)
f (x) d x = [7 f (x0 ) + 32 f (x1 ) + 12 f (x2 ) + 32 f (x3 ) + 7 f (x4 )] − f (ξ ),
x0 45 945
donde x0 < ξ < x4 . (4.28)

Fórmulas de Newton-Cotes abiertas

Las fórmulas de Newton-Cotes abiertas no incluyen los extremos de [a, b] como nodos. Éstas
utilizan los nodos xi = x0 + i h , para cada i = 0, 1,   , n, donde h = (b − a)/(n + 2) y
x0 = a + h. Esto implica que xn = b − h, por lo que etiquetamos los extremos al estable-
cer x−1 = a y xn+1 = bFRPRVHPXHVWUDHQODÀJXUD/DVIyUPXODVDELHUWDVFRQWLHQHQ
todos los nodos que se usan para la aproximación dentro del intervalo abierto (a, b). Las
fórmulas se convertirán en
b xn+1 n
f (x) d x = f (x) d x ≈ ai f (xi ),
a x−1 i=0

b
donde ai = L i (x) d x.
a

Figura 4.6

y = Pn(x)
y = f (x)

a 5 x21 x0 x1 x2 xn xn11 5 b x

El siguiente teorema es análogo al teorema 4.2; su demostración se encuentra en [IK],

p. 314.
n
Teorema 4.3 Suponga que i=0 ai f (xi ) denota la fórmula abierta de (n 1 1) puntos de Newton-Cotes
con x−1 = a, xn+1 = b, y h = (b − a)/(n + 2). Existe ξ ∈ (a, b) para el que
b n n+1
h n+3 f (n+2) (ξ )
f (x) d x = ai f (xi ) + t 2 (t − 1) · · · (t − n) dt,
a i=0
(n + 2)! −1
 4.3 Elementos de integración numérica 149

si n es par y f ∈ C n+2 [a, b], y

b n n+1
h n+2 f (n+1) (ξ )
f (x) d x = ai f (xi ) + t (t − 1) · · · (t − n) dt,
a i=0
(n + 1)! −1

si n es impar y f ∈ C n+1 [a, b].

Observe que, como en el caso de métodos cerrados, tenemos el grado de precisión com-
parativamente superior para los métodos pares que para los métodos impares.
Algunas de las fórmulas de Newton-Cotes abiertas comunes con sus términos de error
son las siguientes:

n 5 0: regla del punto medio

x1
h3
f (x) d x = 2h f (x0 ) + f (ξ ), donde x−1 < ξ < x1 . (4.29)
x−1 3

n 5 1:

x2
3h 3h 3
f (x) d x = [ f (x0 ) + f (x1 )] + f (ξ ), donde x−1 < ξ < x2 . (4.30)
x−1 2 4

n 5 2:

x3
4h 14h 5 (4)
f (x) d x = [2 f (x0 ) − f (x1 ) + 2 f (x2 )] + f (ξ ), (4.31)
x−1 3 45
donde x−1 < ξ < x3 .

n 5 3:

x4
5h 95 5 (4)
f (x) d x = [11 f (x0 ) + f (x1 ) + f (x2 ) + 11 f (x3 )] + h f (ξ ), (4.32)
x−1 24 144
donde x−1 < ξ < x4 .

Ejemplo 2 Compare los resultados de las fórmulas cerradas y abiertas de Newton-Cotes como la ecua-
ción (4.25) a la (4.28) y la ecuación (4.29) a la (4.32) para aproximar
π/4 √
sen x d x = 1 − 2/2 ≈ 0.29289322.
0

Solución Para las fórmulas cerradas, tenemos

(π/4) π
n=1: sen 0 + sen ≈ 0.27768018
2 4
(π/8) π π
n=2: sen 0 + 4 sen + sen ≈ 0.29293264
3 8 4
3(π/12) π π π
n=3: sen 0 + 3 sen + 3 sen + sen ≈ 0.29291070
8 12 6 4
2(π/16) π π 3π π
n=4: 7 sen 0 + 32 sen + 12 sen + 32 sen + 7 sen ≈ 0.29289318
45 16 8 16 4
150 CAPÍTULO 4 Diferenciación numérica e integración

y para las fórmulas abiertas, tenemos

π
n = 0 : 2(π/8) sen ≈ 0.30055887
8
3(π/12) π π
n=1: sen + sen ≈ 0.29798754
2 12 6
4(π/16) π π 3π
n=2: 2 sen − sen + 2 sen ≈ 0.29285866
3 16 8 16
5(π/20) π π 3π π
n=3: 11 sen + sen + sen + 11 sen ≈ 0.29286923
24 20 10 20 5

La tabla 4.8 resume los resultados y muestra los errores de aproximación.

Tabla 4.8 n 0 1 2 3 4
Fórmulas cerradas 0.27768018 0.29293264 0.29291070 0.29289318
Error 0.01521303 0.00003942 0.00001748 0.00000004
Fórmulas abiertas 0.30055887 0.29798754 0.29285866 0.29286923
Error 0.00766565 0.00509432 0.00003456 0.00002399

La sección Conjunto de ejercicios 4.3 está disponible en línea. Encuentre la ruta de

acceso en las páginas preliminares.

4.4 Integración numérica compuesta

En general, el uso de fórmulas de Newton-Cotes es inapropiado sobre largos intervalos de in-
WHJUDFLyQ6HUHTXHULUtDQIyUPXODVGHJUDGRVXSHULRU\ORVYDORUHVGHORVFRHÀFLHQWHVHQHVWDV
fórmulas son difíciles de obtener. Además, las fórmulas de Newton-Cotes están basadas en
los polinomios de interpolación que utilizan nodos igualmente espaciados, un procedimiento
inapropiado sobre intervalos largos debido a la naturaleza oscilatoria de los polinomios de
orden superior.
A menudo, la aproximación por
tramos es efectiva. Recuerde que
En esta sección, analizamos el enfoque por tramos (o fragmentario) para la integración
esto se usó para interpolación de numérica que usa las fórmulas de Newton-Cotes de bajo orden. Estas son las técnicas que se
spline. aplican más a menudo.

Use la regla de Simpson para aproximar 0 e x d x y compare esto con los resultados obteni-
4
Ejemplo 1
dos mediante la suma de las aproximaciones de Simpson para 0 e x dx y 2 e x d x y al sumar
2 4
1 x 2 x 3 x 4 x
éstas con 0 e d x, 1 e d x, 2 e d x, y 3 e d x.

Solución La regla de Simpson en [0, 4] con h 5 2 da

4
2 0
ex d x ≈ (e + 4e2 + e4 ) = 56.76958.
0 3

La respuesta correcta en este caso es e4 − e0 = 53.59815, y el error 23.17143 es mucho

más grande de lo que aceptaríamos normalmente.
 4.4 Integración numérica compuesta 151

Al aplicar la regla de Simpson en cada uno de los intervalos [0, 2] y [2, 4] con h 5 1 da
4 2 4
ex d x = ex d x + ex d x
0 0 2
1 0 1 2
≈ e + 4e + e2 + e + 4e3 + e4
3 3
1 0
= e + 4e + 2e2 + 4e3 + e4
3
= 53.86385.

El error se ha reducido a 20.26570.

Para las integrales en [0, 1], [1, 2], [3, 4], y [3, 4], utilizamos la regla de Simpson cuatro
veces con h = 12, con lo que obtenemos
4 1 2 3 4
ex d x = ex d x + ex d x + ex d x + ex d x
0 0 1 2 3
1 1
≈ e0 + 4e1/2 + e + e + 4e3/2 + e2
6 6
1 2 1 3
+ e + 4e5/2 + e3 + e + 4e7/2 + e4
6 6
1 0
= e + 4e1/2 + 2e + 4e3/2 + 2e2 + 4e5/2 + 2e3 + 4e7/2 + e4
6
= 53.61622.

El error para esta aproximación se ha reducido a 20.01807.

b
Para generalizar este procedimiento para una integral arbitraria a f (x) d x, seleccione
un entero par n. Subdivida el intervalo [a, b] en n subintervalos y aplique la regla de Simpson
en cada par consecutivo de subintervalos (véase ODÀJXUD 

Figura 4.7
y

y 5 f (x)

a 5 x0 x2 x2j22 x2j21 x2j b 5 xn x

Con h = (b − a)/n y x j = a + j h, para cada j = 0, 1, . . . , n, tenemos

b n/2 x2 j
f (x) d x = f (x) d x
a j=1 x2 j−2

n/2
h h 5 (4)
= [ f (x2 j−2 ) + 4 f (x2 j−1 ) + f (x2 j )] − f (ξ j ) ,
j=1
3 90
152 CAPÍTULO 4 Diferenciación numérica e integración

para algunos ξ j con x2 j−2 < ξ j < x2 j, siempre que f ∈ C [a, b]. A través del hecho de
4

que para cada j = 1, 2, . . . , (n/2) − 1 tenemos f (x2 j ) TXH ÀJXUD HQ HO WpUPLQR FRUUHV-
pondiente al intervalo [x2 j−2 , x2 j ] y también en el término correspondiente al intervalo
[x2 j , x2 j+2 ], podemos reducir esta suma a

⎡ ⎤
b (n/2)−1 n/2 n/2
h⎣ h5
f (x) d x = f (x0 ) + 2 f (x2 j ) + 4 f (x2 j−1 ) + f (xn )⎦ − f (4) (ξ j ).
a 3 j=1 j=1
90 j=1

El error relacionado con esta aproximación es

n/2
h5
E( f ) = − f (4) (ξ j ),
90 j=1

donde x2 j−2 < ξ j < x2 j , para cada j = 1, 2, . . . , n/2.

Si f ∈ C 4 [a, b], el teorema de valor extremo 1.9 implica que f (4) asume su máximo y su
mínimo en [a, b]. Puesto que

mín f (4) (x) ≤ f (4) (ξ j ) ≤ máx f (4) (x),

x∈[a,b] x∈[a,b]

tenemos
n/2
n n
mín f (4) (x) ≤ f (4) (ξ j ) ≤ máx f (4) (x)
2 x∈[a,b] j=1
2 x∈[a,b]

y
n/2
2
mín f (4) (x) ≤ f (4) (ξ j ) ≤ máx f (4) (x).
x∈[a,b] n j=1
x∈[a,b]

Por medio del teorema de valor intermedio 1.11, existe m ∈ (a, b) tal que
n/2
2
f (4) (μ) = f (4) (ξ j ).
n j=1

Por lo tanto,
n/2
h5 h5
E( f ) = − f (4) (ξ j ) = − n f (4) (μ),
90 j=1
180

Y, puesto que h = (b − a)/n,

(b − a) 4 (4)
E( f ) = − h f (μ).
180

Estas observaciones producen el siguiente resultado.

Teorema 4.4 Si f ∈ C 4 [a, b], n es par, h = (b − a)/n, y x j = a + j h, para cada j = 0, 1, . . . , n. Exis-

te μ ∈ (a, b) para los que la regla compuesta de Simpson para n subintervalos se puede
reescribir con su término de error como
⎡ ⎤
b (n/2)−1 n/2
h⎣ b − a 4 (4)
f (x) d x = f (a) + 2 f (x2 j ) + 4 f (x2 j−1 ) + f (b)⎦ − h f (μ).
a 3 j=1 j=1
180
 4.4 Integración numérica compuesta 153

Observe que el término de error para la regla compuesta de Simpson es O(h4), mientras
que era O(h5) para la regla estándar de Simpson. Sin embargo, estos índices no son compara-
bles ya que para la regla estándar de Simpson, tenemos h ÀMDHQh = (b − a)/2, pero para la
regla compuesta de Simpson, tenemos h = (b−a)/n, para n un entero par. Esto nos permite
reducir considerablemente el valor de h.
El algoritmo 4.1 utiliza la regla compuesta de Simpson en n subintervalos. Éste es el
algoritmo de cuadratura que se usa con mayor frecuencia para propósito general.

ALGORITMO Regla compuesta de Simpson

4.1 Para aproximar la integral I =
b
f (x) d x:
a

ENTRADA extremos a, b; entero positivo n par.

SALIDA aproximación X I para I.
Paso 1 Haga h = (b − a)/n.
Paso 2 Haga X I 0 = f (a) + f (b);
X I 1 = 0; (Suma de f (x2i−1 ).)
X I 2 = 0. (Suma de f (x2i ).)
Paso 3 Para i = 1, . . . , n − 1 realice los pasos 4 y 5.
Paso 4 Haga X = a + i h.
Paso 5 Si i es par entonces haga X I 2 = X I 2 + f (X )
También determine X I 1 = X I 1 + f (X ).
Paso 6 Haga X I = h(X I 0 + 2 · X I 2 + 4 · X I 1)/3.
Paso 7 SALIDA (X I );
PARE.

El enfoque de subdivisión se puede aplicar a cualquiera de las fórmulas de Newton-

Cotes. Las extensiones de las reglas trapezoidal (véase ODÀJXUD \GHSXQWRPHGLRVHGDQ
sin prueba. La regla trapezoidal sólo requiere un intervalo para cada aplicación, por lo que el
entero n puede ser tanto par como impar.

Figura 4.8

y 5 f (x)

a 5 x0 x1 x j21 xj x n21 b 5 xn x
154 CAPÍTULO 4 Diferenciación numérica e integración

Teorema 4.5 Sean f ∈ C 2 [a, b], h = (b − a)/n, y x j = a + j h, para cada j = 0, 1, . . . , n. Existe m ∈

(a, b) para el que la regla compuesta trapezoidal para n subintervalos se puede reescribir
con este término de error como
⎡ ⎤
b n−1
h⎣ b−a 2
f (x) d x = f (a) + 2 f (x j ) + f (b)⎦ − h f (μ).
a 2 j=1
12

Para la regla compuesta de punto medio, n debe ser, nuevamente, par (véase ODÀJXUD
4.9.)

Figura 4.9
y
y 5 f (x)

a 5 x21 x 0 x1 x2j21 x2j x2j11 xn21 x n b 5 x n11 x

Teorema 4.6 Si f ∈ C 2 [a, b], n es par, h = (b − a)/(n + 2), y x j = a + ( j + 1)h para cada j 5
−1, 0, . . . , n + 1. Existe m ∈ (a, b) para el que la regla compuesta de punto medio para n 1 2
subintervalos se puede reescribir con su término de error como
b n/2
b−a 2
f (x) d x = 2h f (x2 j ) + h f (μ).
a j=0
6

Ejemplo 2 Determine valores de h que garantizarán un error de aproximación menor que 0.00002 al
π
aproximar 0 sen x d x y usar
a) La regla compuesta trapezoidal y b) la regla compuesta de Simpson.

Solución a) La forma de error para la regla trapezoidal compuesta para f (x) = sen x en [0, π ]
es

π h2 π h2 π h2
f (μ) = (− sen μ) = | sen μ|.
12 12 12

3DUDJDUDQWL]DUVXÀFLHQWHSUHFLVLyQFRQHVWDWpFQLFDQHFHVLWDPRVWHQHU

π h2 π h2
| sen μ| ≤ < 0.00002.
12 12

Puesto que h = π/n, necesitamos

1/2
π3 π3
< 0.00002, lo cual implica que n > ≈ 359.44,
12n 2 12(0.00002)
 4.4 Integración numérica compuesta 155

y la regla compuesta trapezoidal requiere n ≥ 360.

b) La forma de error para la regla compuesta de Simpson para f (x) = sen x en [0, π ] es

π h 4 (4) π h4 π h4
f (μ) = sen μ = | sen μ|.
180 180 180

3DUDJDUDQWL]DUVXÀFLHQWHSUHFLVLyQFRQHVWDWpFQLFDQHFHVLWDPRVWHQHU

π h4 π h4
| sen μ| ≤ < 0.00002.
180 180

A través del hecho de que n = π/ h nos da

1/4
π5 π5
< 0.00002, lo cual implica que n > ≈ 17.07.
180n 4 180(0.00002)

Por lo tanto, la regla compuesta de Simpson sólo requiere n ≥ 18.

La regla compuesta de Simpson con n 5 18 nos da
⎡ ⎤
π 8 9
π ⎣ jπ (2 j − 1)π ⎦
sen x d x ≈ 2 sen +4 sen = 2.0000104.
0 54 j=1
9 j=1
18

Esto es preciso dentro de aproximadamente 1025 porque el valor verdadero es − cos(π) −

(− cos(0)) = 2.
La regla compuesta de Simpson es la selección clara si desea minimizar los cálculos.
Para propósitos de comparación, considere la regla compuesta trapezoidal por medio de
h = π/18 para la integral en el ejemplo 2. Esta aproximación utiliza las mismas evalua-
ciones de función que la regla compuesta de Simpson, pero la aproximación en este caso,
⎡ ⎤ ⎡ ⎤
π 17 17
π ⎣ jπ π jπ
sen x d x ≈ 2 sen + sen 0 + sen π ⎦ = ⎣2 sen ⎦
0 36 j=1
18 36 j=1
18

= 1.9949205,

sólo es precisa para aproximadamente 5 3 1023.

Estabilidad del error de redondeo

En el ejemplo 2, observamos que garantizar una precisión de 2 3 1025 para aproximar
π
0 sen x d x requería 360 subdivisiones de [0, /] para la regla compuesta trapezoidal y sólo 18
para la regla compuesta de Simpson. Además del hecho de que se necesitan menos cálcu-
los para la técnica de Simpson, usted podría sospechar que este método también implica
menos error de redondeo. Sin embargo, una propiedad importante compartida por todas
las técnicas de integración compuesta es una estabilidad respecto al error de redondeo.
Es decir, el error de redondeo no depende del número cálculos realizados.
Se espera que la integración
numérica sea estable, mientras
Para demostrar este hecho considerablemente sorprendente, suponga que aplicamos la
que la diferenciación numérica es regla compuesta de Simpson con n subintervalos a una función f en [a, b] y determine
inestable. la máxima cota para el error de redondeo. Suponga que f (xi ) se aproxima mediante f˜(xi )
y que

f (xi ) = f˜(xi ) + ei , para cada i = 0, 1, . . . , n,

156 CAPÍTULO 4 Diferenciación numérica e integración

donde ei denota el error de redondeo asociado con el uso de f˜(xi ) para aproximar f (xi ).
Entonces, el error acumulado, e(h), en la regla compuesta de Simpson es

⎡ ⎤
(n/2)−1 n/2
h⎣
e(h) = e0 + 2 e2 j + 4 e2 j−1 + en ⎦
3 j=1 j=1
⎡ ⎤
(n/2)−1 n/2
h⎣
≤ |e0 | + 2 |e2 j | + 4 |e2 j−1 | + |en |⎦ .
3 j=1 j=1

Si los errores de redondeo están limitados de manera uniforme por ε, entonces

h n n h
e(h) ≤ ε+2 −1 ε+4 ε + ε = 3nε = nhε.
3 2 2 3

Sin embargo, nh = b − a , por lo que

e(h) ≤ (b − a)ε,

una cota independiente de h (y n (VWRVLJQLÀFDTXHDXQTXHTXL]iQHFHVLWHPRVGLYLGLUXQ

intervalo en más partes para garantizar precisión, el cálculo incrementado que se requiere
no aumenta el error de redondeo. Este resultado implica que el procedimiento es estable
conforme h se aproxima a cero. Recuerde que esto no es verdad para los procedimientos de
diferenciación numérica considerados al principio de este capítulo.

La sección Conjunto de ejercicios 4.4 está disponible en línea. Encuentre la ruta de

acceso en las páginas preliminares.

4.5 Integración de Romberg

En esta sección ilustraremos la forma en la que la extrapolación de Richardson aplicada a los
resultados de la regla compuesta trapezoidal se puede usar para obtener aproximaciones de
alta precisión con poco costo computacional.
En la sección 4.4 encontramos que la regla compuesta trapezoidal tiene un error de trun-
camiento de orden O(h2 (QHVSHFtÀFRPRVWUDPRVTXHSDUDh = (b − a)/n y x j = a + j h,
tenemos
⎡ ⎤
b n−1
h⎣ (b − a) f (μ) 2
f (x) d x = f (a) + 2 f (x j ) + f (b)⎦ − h ,
a 2 j=1
12

para algún numero m en (a, b).

Con un método alternativo, se puede mostrar (consulte [RR], pp. 136–140) que si
f ∈ C ∞ [a, b], la regla compuesta trapezoidal también se puede escribir con un término de
error de la forma
⎡ ⎤
b n−1
h⎣
f (x) d x = f (a) + 2 f (x j ) + f (b)⎦ + K 1 h 2 + K 2 h 4 + K 3 h 6 + · · · ,
a 2 j=1 (4.33)

donde Ki es una constante que depende solamente de f (2i−1) (a) y f (2i−1) (b).
 4.5 Integración de Romberg 157

Recuerde, de la sección 4.2, que la extrapolación de Richardson se puede realizar en

cualquier procedimiento de aproximación cuyo error de truncamiento es de la forma
m−1
K j h α j + O(h αm ),
j=1

para un conjunto de constantes Kj y donde α1 < α2 < α3 < · · · < αm . En esa sección reali-
zamos demostraciones para ilustrar qué tan efectiva es esta técnica cuando el procedimiento
de aproximación tiene un error de truncamiento sólo con potencias pares de h, es decir,
cuando el error de truncamiento tiene la forma
m−1
K j h 2 j + O(h 2m ).
j=1

Werner Romberg (1909–2093) Puesto que la regla trapezoidal compuesta tiene esta forma, es un candidato obvio para extra-
concibió este procedimiento para polación. Esto resulta en una técnica conocida como integración de Romberg.
b
mejorar la precisión de la regla Para aproximar la integral a f (x) d x, utilizamos los resultados de la regla compues-
trapezoidal al eliminar términos ta trapezoidal con n = 1, 2, 4, 8, 16, . . . , e indicamos las aproximaciones resultantes, res-
sucesivos en la expansión
asintótica en 1955. pectivamente mediante R1,1, R2,1, R3,1, y así sucesivamente. A continuación, aplicamos la
extrapolación de la forma determinada en la sección 4.2; es decir, obtenemos O(h4) aproxi-
maciones R2,2, R3,2, R4,2, y así sucesivamente, por medio de

1
Rk,2 = Rk,1 + (Rk,1 − Rk−1,1 ), para k = 2, 3, . . .
3

y, entonces, las O(h6) aproximaciones R3,3, R4,3, R5,3, se pueden obtener mediante

1
Rk,3 = Rk,2 + (Rk,2 − Rk−1,2 ), para k = 3, 4, . . .
15

En general, después de que se han obtenido las aproximaciones Rk, j−1 adecuadas, determi-
namos las aproximaciones O(h 2 j ) a partir de
1
Rk, j = Rk, j−1 + (Rk, j−1 − Rk−1, j−1 ), para k = j, j + 1, . . .
4 j−1 − 1
π
Ejemplo 1 Use la regla compuesta trapezoidal para encontrar aproximaciones para 0 sen x d x con
n 5 1, 2, 4, 8 y 16. A continuación, realice la integración de Romberg en los resultados.
La regla compuesta trapezoidal para los diferentes valores de n proporciona las siguientes
aproximaciones para el valor verdadero 2:
π
R1,1 = [sen 0 + sen π ] = 0,
2
π π
R2,1 = sen 0 + 2 sen + sen π = 1.57079633,
4 2
π π π 3π
R3,1 = sen 0 + 2 sen + sen + sen + sen π = 1.89611890,
8 4 2 4
π π π 3π 7π
R4,1 = sen 0 + 2 sen + sen + · · · + sen + sen + sen π
16 8 4 4 8
=1.97423160, y
π π π 7π 15π
R5,1 = sen 0 + 2 sen + sen + · · · + sen + sen + sen π
32 16 8 8 16
= 1.99357034.
158 CAPÍTULO 4 Diferenciación numérica e integración

Las aproximaciones O(h4) son

1
R2,2 = R2,1 + (R2,1 − R1,1 ) = 2.09439511,
3
1
R3,2 = R3,1 + (R3,1 − R2,1 ) = 2.00455976,
3
1
R4,2 = R4,1 + (R4,1 − R3,1 ) = 2.00026917, y
3
1
R5,2 = R5,1 + (R5,1 − R4,1 ) = 2.00001659.
3

Las aproximaciones O(h6) son

1
R3,3 = R3,2 + (R3,2 − R2,2 ) = 1.99857073,
15
1
R4,3 = R4,2 + (R4,2 − R3,2 ) = 1.99998313, y
15
1
R5,3 = R5,2 + (R5,2 − R4,2 ) = 1.99999975.
15

Las dos aproximaciones O(h8) son

1 1
R4,4 = R4,3 + (R4,3 − R3,3 ) = 2.00000555 y R5,4 = R5,3 + (R5,3 − R4,3 )
63 63
= 2.00000001,

y la aproximación O(h10 ÀQDOHV

1
R5,5 = R5,4 + (R5,4 − R4,4 ) = 1.99999999.
255
Estos resultados se muestran en la tabla 4.9.
Tabla 4.9 0
1.57079633 2.09439511
1.89611890 2.00455976 1.99857073
1.97423160 2.00026917 1.99998313 2.00000555
1.99357034 2.00001659 1.99999975 2.00000001 1.99999999

Observe que al generar las aproximaciones para la regla compuesta trapezoidal en el

ejemplo 1, cada aproximación consecutiva incluía todas las evaluaciones de la función des-
de la aproximación previa. Es decir, R1,1 utilizaba evaluaciones en 0 y π, y R2,1 usaba estas
evaluaciones y añadía una evaluación en el punto intermedio π/2. Entonces R3,1 utilizaba
las evaluaciones de R2,1 y añadía dos intermedios adicionales en π/4 y 3π/4. Este patrón
continúa con R4,1 mediante las mismas evaluaciones como R3,1, pero al añadir evaluaciones
en los cuatro puntos intermedios π/8, 3π/8, 5π/8, 7π/8, y así sucesivamente.
Este procedimiento de evaluación para las aproximaciones de la regla compuesta tra-
dicional se mantiene para una integral en cualquier intervalo [a, b]. En general, la regla
compuesta trapezoidal denotaba Rk+1,1 utiliza las mismas evaluaciones que Rk,1, pero añade
evaluaciones en los puntos intermedios 2 k−23RUORWDQWRHOFiOFXORHÀFLHQWHGHHVWDVDSUR-
ximaciones puede realizarse de manera recursiva.
b
Para obtener las aproximaciones de la regla compuesta trapezoidal para a f (x) d x,
k−1
entonces h k = (b − a)/m k = (b − a)/2 .

h1 (b − a)
R1,1 = [ f (a) + f (b)] = [ f (a) + f (b)],
2 2
 4.5 Integración de Romberg 159

h2
R2,1 = [ f (a) + f (b) + 2 f (a + h 2 )].
2
Al reescribir este resultado para R2,1, podemos incorporar la aproximación previamente
determinada R1,1

(b − a) (b − a) 1
R2,1 = f (a) + f (b) + 2 f a+ = [R1,1 + h 1 f (a + h 2 )].
4 2 2

De manera similar, podemos escribir

1
R3,1 = {R2,1 + h 2 [ f (a + h 3 ) + f (a + 3h 3 )]},
2

y, en general (véaseODÀJXUD WHQHPRV

⎡ ⎤
2k−2
1
Rk,1 = ⎣ Rk−1,1 + h k−1 f (a + (2i − 1)h k )⎦ , (4.34)
2 i=1

para cada k = 2, 3, . . . , n. (Consulte los ejercicios 18 y 19.)

Figura 4.10

y y y

y 5 f (x) y 5 f (x) y 5 f (x)

R1,1 R 2,1 R 3,1

a b x a b x a b x

Entonces, la extrapolación se usa para producir aproximaciones O(h 2k j ) mediante

1
Rk, j = Rk, j−1 + (Rk, j−1 − Rk−1, j−1 ), para k = j, j + 1, . . . ,
4 j−1 −1
como se muestra en la tabla 4.10.

Tabla 4.10 k O h 2k O h 4k O h 6k O h 8k O h 2n
k

1 R1,1
2 R2,1 R2,2
3 R3,1 R3,2 R3,3
4 R4,1 R4,2 R4,3 R4,4
.. .. .. .. .. ..
. . . . . .
n Rn,1 Rn,2 Rn,3 Rn,4 ··· Rn,n
160 CAPÍTULO 4 Diferenciación numérica e integración

El método efectivo para construir la tabla de Romberg utiliza el orden más alto de la
DSUR[LPDFLyQHQFDGDSDVR(VGHFLUFDOFXODODVHQWUDGDVÀODSRUÀODHQRUGHQR1,1, R2,1, R2,2,
R3,1, R3,2, R3,3, y así sucesivamente. Esto también permite calcular unaQXHYDÀODFRPSOHWD
en la tabla al realizar sólo una aplicación adicional de la regla compuesta trapezoidal. A con-
tinuación, usa un promediado simple de los valores previamente calculados para obtener las
HQWUDGDVUHVWDQWHVHQODÀOD5HFXHUGH

‡ &DOFXODUODWDEODGH5RPEHUJXQDÀODFRPSOHWDDODYH]

π
Ejemplo 2 $xDGDXQDÀODGHH[WUDSRODFLyQDGLFLRQDODODWDEODSDUD 0 sen x d x.

Solución 3DUDREWHQHUODÀODDGLFLRQDOQHFHVLWDPRVODDSUR[LPDFLyQWUDSH]RLGDO

⎡ ⎤
24
1⎣ π (2k − 1)π ⎦
R6,1 = R5,1 + sen = 1.99839336.
2 16 k=1
32

El valor en la tabla 4.9 da

1 1
R6,2 = R6,1 + (R6,1 − R5,1 ) = 1.99839336 + (1.99839336 − 1.99357035)
3 3
= 2.00000103,
1 1
R6,3 = R6,2 + (R6,2 − R5,2 ) = 2.00000103 + (2.00000103 − 2.00001659)
15 15
= 2.00000000,
1 1
R6,4 = R6,3 + (R6,3 − R5,3 ) = 2.00000000, R6,5 = R6,4 + (R6,4 − R5,4 )
63 255
= 2.00000000,

y R6,6 = R6,5 + 10231

(R6,5 − R5,5 ) = 2.00000000. La nueva tabla de extrapolación se mues-
tra en la tabla 4.11.

Tabla 4.11 0
1.57079633 2.09439511
1.89611890 2.00455976 1.99857073
1.97423160 2.00026917 1.99998313 2.00000555
1.99357034 2.00001659 1.99999975 2.00000001 1.99999999
1.99839336 2.00000103 2.00000000 2.00000000 2.00000000 2.00000000

2EVHUYHTXHWRGRVORVYDORUHVH[WUDSRODGRVH[FHSWRSRUHOSULPHUR HQODSULPHUDÀODGH
la segunda columna) son más precisos que la mejor aproximación compuesta trapezoidal (en la
~OWLPDÀODGHODSULPHUDFROXPQD $SHVDUGHTXHH[LVWHQHQWUDGDVHQODWDEODVólo
la sexta en la columna izquierda requiere evaluaciones de función ya que son las únicas en-
tradas generadas por la regla compuesta trapezoidal; las otras entradas se obtienen mediante
un proceso de promediado. De hecho, debido a la relación de recurrencia de los términos
en la columna izquierda, las únicas evaluaciones de función son aquellas para calcular la
aproximación de la regla compuesta trapezoidal. En general, Rk,1 requiere evaluaciones de
función 1 + 2k−1, por lo que en este caso se necesita 1 + 25 = 33.
 4.5 Integración de Romberg 161

El algoritmo 4.2 usa el procedimiento recursivo para encontrar las aproximaciones de la

UHJODFRPSXHVWDWUDSH]RLGDO\FDOFXODORVUHVXOWDGRVHQODWDEODÀODSRUÀOD

ALGORITMO Integración de Romberg

b
4.2 Para aproximar la integral I = f (x) d x, seleccione un entero n > 0.
a

ENTRADA extremos a, b; entero n.

SALIDA un arreglo R. (Calcule R por filas; sólo se almacenan las últimas dos filas.)
Paso 1 Haga h = b − a;
R1,1 = h2 ( f (a) + f (b)).
Paso 2 SALIDA (R1,1 ).
Paso 3 Para i = 2, . . . , n realice los pasos 4–8.
⎡ ⎤
2i−2
1⎣
Paso 4 Haga R2,1 = R1,1 + h f (a + (k − 0.5)h)⎦.
2 k=1

(Aproximación a partir del método trapezoidal.)

Paso 5 Para j = 2, . . . , i
R2, j−1 − R1, j−1
haga R2, j = R2, j−1 + . (Extrapolación.)
4 j−1 − 1
Paso 6 SALIDA (R2, j para j = 1, 2, . . . , i).
Paso 7 Haga h = h/2.
Paso 8 Para j = 1, 2, . . . , i determine R1, j = R2, j . (Actualice la fila 1 de R.)
Paso 9 PARE.

El algoritmo 4.2 requiere un entero preestablecido n para determinar el número de

ÀODV TXH VH YD D JHQHUDU7DPELpQ HVWDEOHFHUHPRV XQD WROHUDQFLD GH HUURU SDUD OD DSUR[L-
mación y generar n, dentro de algúna cota superior hasta que las entradas diagonales con-
secutivas Rn−1,n−1 y Rn,n concuerden dentro de la tolerancia. Para evitar la posibilidad de
TXHGRVHOHPHQWRVGHÀODFRQVHFXWLYRVFRQFXHUGHQXQRVFRQRWURVSHURQRFRQHOYDORUGH
la integral que se está aproximando, es común generar aproximaciones hasta que no sólo
|Rn−1,n−1 − Rn,n | esté dentro de la tolerancia, sino también |Rn−2,n−2 − Rn−1,n−1 |. A pesar
de que no es una protección universal, esto garantizará que dos conjuntos de aproximaciones
JHQHUDGRVGHPDQHUDGLIHUHQWHFRQFXHUGHQGHQWURGHODWROHUDQFLDHVSHFLÀFDGDDQWHVGHRn,n
VHDFHSWHFRPRVXÀFLHQWHPHQWHSUHFLVR
La integración de Romberg aplicada a una función f en el intervalo [a, b] depende de la
suposición de que la regla compuesta trapezoidal tienen un término de error que se puede
expresar en la forma de la ecuación (4.33); es decir, debemos tener f ∈ C 2k+2 [a, b] para la
El adjetivo cauteloso que se usa k-ésimaÀODTXHVHYDDJHQHUDU/RVDOJRULWPRVGHSURSyVLWRJHQHUDOSRUPHGLRGHODLQWH-
en la descripción de un método JUDFLyQGH5RPEHUJLQFOX\HQXQDYHULÀFDFLyQHQFDGDHWDSDSDUDJDUDQWL]DUHOFXPSOLPLHQWR
numérico indica que se incluye de esta suposición. Estos métodos se conocen como algoritmos cautelosos de Romberg y
XQDYHULÀFDFLyQSDUDGHWHUPLQDU
si las hipótesis de continuidad
se describen en [Joh]. Esta referencia también describe métodos para utilizar la técnica de
tienen alguna probabilidad de ser Romberg como un procedimiento adaptable, similar a la regla adaptable de Simpson, la cual
verdaderas. se analizará en la sección 4.6.

La sección Conjunto de ejercicios 4.5 está disponible en línea. Encuentre la ruta de

acceso en páginas preliminares.
162 CAPÍTULO 4 Diferenciación numérica e integración

4.6 Métodos de cuadratura adaptable

Las fórmulas compuestas son muy efectivas en muchas situaciones, pero ocasionalmente su-
fren porque requieren el uso de nodos igualmente espaciados. Esto es inadecuado al integrar
una función sobre un intervalo que contiene ambas regiones con gran variación funcional y
regiones con variación funcional pequeña.

Ilustración La única solución a la ecuación diferencial y + 6y + 25 = 0 que adicionalmente satis-

face y(0) = 0 y y (0) = 4 es y(x) = e−3x sen 4x. Las funciones de este tipo son comunes
en ingeniería mecánica porque describen ciertas características de sistemas absorbentes de
muelle y amortiguador y en ingeniería eléctrica porque son soluciones comunes a los pro-
EOHPDVIXQGDPHQWDOHVGHFLUFXLWRV/DJUiÀFDGHy(x) para x en el intervalo [0, 4] se muestra
HQODÀJXUD

Figura 4.11

0.5

0.4

0.3
y (x) = e23xsen 4x
0.2

0.1
1
2
2 3 4 x
20.1

Suponga que necesitamos la integral de y(x HQ>@/DJUiÀFDLQGLFDTXHODLQWHJUDO

en [3, 4] debe estar muy cerca de 0 y en [2, 3] tampoco se esperaría que fuera grande. Sin
HPEDUJRHQ>@H[LVWHYDULDFLyQVLJQLÀFDWLYDGHODIXQFLyQ\QRHVGHOWRGRFODUDFXiOHVOD
integral en este intervalo. Este es un ejemplo de una situación en donde la integración com-
puesta sería inadecuada. Se podría usar un método de orden bajo en [2, 4], pero se necesitaría
un método de orden superior en [0, 2].

La pregunta que consideraríamos en esta sección es:

• ¿Cómo podemos determinar la técnica que deberíamos aplicar en las diferentes partes del
LQWHUYDORGHLQWHJUDFLyQ\TXpWDQSUHFLVDSRGHPRVHVSHUDUTXHVHDODDSUR[LPDFLyQÀQDO"

Veremos que en ciertas condiciones razonables, podemos responder esta pregunta y también
determinar aproximaciones que satisfacen requisitos de precisión determinados.
Si el error de aproximación para una integral en un intervalo determinado está distri-
buido de manera equitativa, se necesita un tamaño de paso más pequeño para las grandes
UHJLRQHVGHYDULDFLyQTXHSDUDDTXHOODVFRQPHQRVYDULDFLyQ8QDWpFQLFDHÀFLHQWHSDUDHVWH
 4.6 Métodos de cuadratura adaptable 163

tipo de problema debería predecir la cantidad de variación funcional y adaptar el tamaño de

paso conforme sea necesario. Estos métodos reciben el nombre de métodos de cuadratura
adaptable. Los métodos adaptables son especialmente populares para la inclusión en pa-
TXHWHVSURIHVLRQDOHVGHVRIWZDUHSRUTXHDGHPiVGHVHUHÀFLHQWHVHQJHQHUDOSURSRUFLRQDQ
DSUR[LPDFLRQHVTXHVHHQFXHQWUDQGHQWURGHXQDWROHUDQFLDHVSHFtÀFDGHWHUPLQDGD
En esta sección, consideramos un método de cuadratura y veremos cómo se puede usar
para reducir el error de aproximación y también predecir un error calculado para la apro-
ximación que no depende del conocimiento de las derivadas superiores de la función. El
método que analizamos está basado en la regla compuesta de Simpson, pero la técnica se
PRGLÀFDIiFLOPHQWHSDUDXWLOL]DURWURVSURFHGLPLHQWRVFRPSXHVWRV
b
Suponga que queremos aproximar a f (x) d xGHQWURGHXQDWROHUDQFLDHVSHFtÀFDε > 0.
El primer paso es aplicar la regla de Simpson con tamaño de paso h = (b − a)/2. Esto pro-
duce (véaseODÀJXUD
b
h 5 (4)
f (x) d x = S(a, b) − f (ξ ), para algunos ξ en(a, b), (4.35)
a 90

donde denotamos la aproximación de la regla de Simpson en [a, b] mediante

h
S(a, b) = [ f (a) + 4 f (a + h) + f (b)].
3

Figura 4.12
y

S(a, b) y 5 f (x)

a b x
h h

El siguiente paso es determinar una aproximación de precisión que no requiere f (4) (ξ ).

Para hacerlo aplicamos la regla compuesta de Simpson con n 5 4 y el tamaño de paso
(b−a)/4 = h/2, lo cual nos da
b
h h 3h
f (x) d x = f (a) + 4 f a+ + 2 f (a + h) + 4 f a+ + f (b)
a 6 2 2
4
h (b − a) (4)
− f (ξ̃ ), (4.36)
2 180

para algunos j en (a, b 3DUDVLPSOLÀFDUODQRWDFLyQVL

a+b h h
S a, = f (a) + 4 f a+ + f (a + h)
2 6 2
164 CAPÍTULO 4 Diferenciación numérica e integración

a+b h 3h
S ,b = f (a + h) + 4 f a+ + f (b) .
2 6 2

Entonces, la ecuación (4.36) se puede reescribir (véase la figura 4.13) como

b
a+b a+b 1 h5
f (x) d x = S a, +S ,b − f (4) (ξ̃ ). (4.37)
a 2 2 16 90

Figura 4.13

(
S a,
a1b
2 (
1S
a1b
2
(, b (
y 5 f (x)

a a1b b x
2
h
2

El cálculo de error se deriva al suponer que ξ ≈ ξ̃ o, con mayor precisión, que f (4) (ξ ) ≈
(4)
f (ξ̃ ), y el éxito de la técnica depende de la precisión de esta suposición. Si es precisa,
entonces, al equiparar las integrales en las ecuaciones (4.35) y (4.37) obtenemos

a+b a+b 1 h5 h 5 (4)

S a, +S ,b − f (4) (ξ ) ≈ S(a, b) − f (ξ ),
2 2 16 90 90

por lo que

h 5 (4) 16 a+b a+b

f (ξ ) ≈ S(a, b) − S a, −S ,b .
90 15 2 2

Al utilizar esta técnica en la ecuación (4.37) produce el cálculo de error

b
a+b a+b
f (x) d x − S a, −S ,b
a 2 2

1 h5 1 a+b a+b
≈ f (4) (ξ ) ≈ S(a, b) − S a, −S ,b .
16 90 15 2 2
 4.6 Métodos de cuadratura adaptable 165

b
Esto implica que S(a, (a + b)/2) + S((a + b)/2, b) se aproxima a a f (x) d x aproxima-
damente 15 veces mejor que lo que concuerda con el valor calculado S(a, b). Por lo tanto, si

a+b a+b
S(a, b) − S a, −S ,b < 15ε, (4.38)
2 2

esperamos tener
b
a+b a+b
f (x) d x − S a, −S ,b < ε, (4.39)
a 2 2

a+b a+b
S a, +S ,b
2 2
b
VHVXSRQHTXHHVXQDDSUR[LPDFLyQVXÀFLHQWHPHQWHSUHFLVDSDUD a f (x) d x.

Ejemplo 1 9HULÀTXHODSUHFLVLyQGHOFiOFXORGHHUURUGHWHUPLQDGRSRUODVGHVLJXDOGDGHV  \  

al aplicar la integral

π/2
sen x d x = 1
0

al comparar

π/2
1 π π π π π π π
S 0, − S 0, −S , con sen x d x − S 0, −S , .
15 2 4 4 2 0 4 4 2

Solución 7HQHPRV

π π/4 π π π √
S 0, = sen 0 + 4 sen + sen = (2 2 + 1) = 1.002279878
2 3 4 2 12

π π π π/8 π π 3π π
S 0, +S , = sen 0 + 4 sen + 2 sen + 4 sen + sen
4 4 2 3 8 4 8 2
= 1.000134585.

Por lo que,

π π π π
S 0, − S 0, −S , = |1.002279878 − 1.000134585| = 0.002145293.
2 4 4 2

El cálculo para el error obtenido al utilizar S(a, (a + b)) + S((a + b), b) para aproximar
b
a f (x) d x es, por consiguiente,

1 π π π π
S 0, − S 0, −S , = 0.000143020,
15 2 4 4 2
166 CAPÍTULO 4 Diferenciación numérica e integración

el cual se aproxima de cerca al error real

π/2
sen x d x − 1.000134585 = 0.000134585,
0

aunque Dx4 sen x = sen xYDUtDVLJQLÀFDWLYDPHQWHHQHOLQWHUYDOR(0, π/2).

&XDQGR ODV DSUR[LPDFLRQHV HQ OD GHVLJXDOGDG   GLÀHUHQ SRU PiV GH 15ε, pode-
mos aplicar la técnica de la regla de Simpson de manera individual a los subintervalos
[ a, (a + b)/2] y [(a + b)/2, b]. A continuación, usamos el procedimiento de cálculo de
error para determinar si la aproximación para la integral en cada subintervalo se encuentra
dentro de una tolerancia de ε/2. En este caso, sumamos las aproximaciones para producir
b
una aproximación de a f (x) d x dentro de la tolerancia ε.
Si la aproximación en uno de los subintervalos no se encuentra dentro de la tolerancia
ε/2, entonces ese subintervalo se subdivide en sí mismo y el procedimiento se reaplica a los
dos subintervalos para determinar si la aproximación en cada subintervalo es precisa dentro de
ε/4. Este procedimiento de reducción a la mitad continúa hasta que cada parte está dentro
de la tolerancia requerida.
Se pueden construir problemas para los que esta tolerancia nunca se cumplirá, pero
normalmente la técnica es exitosa porque cada subdivisión aumenta la precisión de la apro-
ximación en un factor de 16 mientras requiere un factor de precisión aumentado en sólo 2.
El algoritmo 4.3 describe detalladamente este procedimiento de cuadratura adaptable
SDUDODUHJODGH6LPSVRQDSHVDUGHTXHODLPSOHPHQWDFLyQGLÀHUHOLJHUDPHQWHGHODQiOLVLV
anterior. Por ejemplo, en el paso 1, la tolerancia se ha establecido en 10ε en lugar de en 15ε en
la desigualdad (4.38). Esta cota se elige de manera conservadora para compensar el error
en la suposición f (4) (ξ ) ≈ f (4) (ξ̃ ). En problemas donde se sabe que f (4) varía ampliamente,
esta cota debería disminuir todavía más.
Es buena idea incluir un margen El procedimiento listado en el algoritmo, primero aproxima la integral en el subinter-
de seguridad cuando es imposible YDORVLWXDGRPiVDODL]TXLHUGDHQXQDVXEGLYLVLyQ(VWRUHTXLHUHDOPDFHQDPLHQWRHÀFLHQWH
YHULÀFDUODVVXSRVLFLRQHVGH
precisión. y recordar las evaluaciones funcionales que se han calculado antes para los nodos en los
subintervalos situados a la mitad derecha. Los pasos 3, 4 y 5 contienen un procedimiento de
apilamiento con un indicador para seguir los datos requeridos para calcular la aproximación
en el subintervalo inmediatamente adyacente y a la derecha del subintervalo en el que se
genera la aproximación. El método es más fácil de implementar por medio de lenguaje de
programación recursivo.

ALGORITMO Cuadratura adaptable

4.3 b
Para aproximar la integral I = f (x) d x dentro de una tolerancia determinada:
a

ENTRADA extremos a, b; tolerancia TOL; cota N para diferentes niveles.

SALIDA aproximación APP o mensaje N superado.
Paso 1 Haga APP = 0;
i = 1;
TOLi = 10 TOL;
ai = a;
h i = (b − a)/2;
FAi = f (a);
FCi = f (a + h i );
FBi = f (b);
Si = h i (FAi + 4FCi + FBi )/3; (Aproximación a partir del método
de Simpson para intervalo completo.)
L i = 1.
Paso 2 Si i > 0 haga los pasos 3–5.
 4.6 Métodos de cuadratura adaptable 167

Paso 3 Haga FD = f (ai + h i /2);

FE = f (ai + 3h i /2);
S1 = h i (FAi + 4FD + FCi )/6; (Aproximaciones para el método
de Simpson para mitades de
subintervalos.)
S2 = h i (FCi + 4FE + FBi )/6;
v1 = ai ; (Guardar datos en este nivel.)
v2 = FAi ;
v3 = FCi ;
v4 = FBi ;
v5 = h i ;
v6 = TOLi ;
v7 = Si ;
v8 = L i .
Paso 4 Haga i = i − 1. (Borrar el nivel.)
Paso 5 Si|S1 + S2 − v7 | < v6
entonces determine APP = APP + (S1 + S2)
también
si (v8 ≥ N )
entonces
SALIDA ('NIVEL EXCEDIDO'); (Falla del procedimiento.)
PARE.
también (Añadir un nivel. )
haga i = i + 1; (Datos para subintervalo mitad derecha.)
ai = v1 + v5 ;
FAi = v3 ;
FCi = FE;
FBi = v4 ;
h i = v5 /2;
TOLi = v6 /2;
Si = S2;
L i = v8 + 1;
haga i = i + 1; (Datos para subintervalo mitad izquierda.)
ai = v1 ;
FAi = v2 ;
FCi = FD;
FBi = v3 ;
h i = h i−1 ;
TOLi = TOLi−1 ;
Si = S1;
L i = L i−1 .
Paso 6 SALIDA (APP); (APP se aproxima a I dentro de TOL. )
PARE.

Ilustración /DJUiÀFDGHODIXQFLyQ f (x) = (100/x 2 ) sen(10/x) para xHQ>@VHPXHVWUDHQODÀ-

gura 4.13. Por medio del algoritmo de cuadratura adaptable 4.3 con tolerancia 1024 para
aproximar 1 f (x) d x produce 21.426014, un resultado preciso dentro de 1.1 × 10−5.
3

La aproximación requería que la regla de Simpson con n 5 4 se realice sobre los 23 subin-
WHUYDORVFX\RVH[WUHPRVVHPXHVWUDQHQHOHMHKRUL]RQWDOHQODÀJXUD(OQ~PHURWRWDOGH
evaluaciones funcionales requerido para esta aproximación es 93.
168 CAPÍTULO 4 Diferenciación numérica e integración

Figura 4.14

60

50

40

30 y = f (x) =
100
x2
sen( (
10
x
20

10
2.75 3.0

1.0 1.25 1.5 1.75 2.0 2.25 2.5 x

210

220

230

240

250

260

El valor más grande de h para el que la regla compuesta de Simpson estándar provee
precisión dentro de 1024 es h = 1/88. Esta aplicación requiere 177 evaluaciones de función,
aproximadamente dos veces más que la cuadratura adaptable.

La sección Conjunto de ejercicios 4.6 está disponible en línea: Encuentre la ruta de

acceso en las páginas preliminares.

4.7 Cuadratura gaussiana

Las fórmulas de Newton-Cotes en la sección 4.3 se dedujeron al integrar polinomios de in-
terpolación. El término de error en el polinomio interpolante de grado n implica la derivada
(n 1 1) de la función que se va a aproximar, por lo que la fórmula Newton-Cotes es exacta
al aproximar la integral de cualquier polinomio de grado menor o igual a n.
7RGDVODVIyUPXODVGH1HZWRQ&RWHVXWLOL]DQYDORUHVGHODIXQFLyQHQSXQWRVLJXDOPHQWH
espaciados. Esta restricción es conveniente cuando se combinan las fórmulas para formar las
UHJODVFRPSXHVWDVTXHFRQVLGHUDPRVHQODVHFFLyQSHURSXHGHGLVPLQXLUVLJQLÀFDWLYD-
mente la precisión de la aproximación. Considere, por ejemplo, la regla trapezoidal aplicada
SDUDGHWHUPLQDUODVLQWHJUDOHVGHODVIXQFLRQHVFX\DVJUiÀFDVVHPXHVWUDQHQODÀJXUD
 4.7 Cuadratura gaussiana 169

Figura 4.15

y y y

y 5 f (x) y 5 f (x)
y 5 f (x)

a 5 x1 x2 5 b x a 5 x1 x2 5 b x a 5 x1 x2 5 b x

La regla trapezoidal aproxima la integral de la función al integrar la función lineal que

XQHORVH[WUHPRVGHODJUiÀFDGHODIXQFLyQ3HURSUREDEOHPHQWHQRHVODPHMRUOtQHDSDUD
DSUR[LPDUODLQWHJUDO(QPXFKRVFDVRVODVOtQHDVFRPRODVPRVWUDGDVHQODÀJXUDSUR-
bablemente proporcionarían mucho mejores aproximaciones.

Figura 4.16

y y y

y 5 f (x)
y 5 f (x)
y 5 f (x)

a x1 x2 b x a x1 x2 b x a x1 x2 b x

En la cuadratura Gaussiana, los puntos para evaluación se seleccionan de manera óptima

en lugar de nodos igualmente espaciados. Los nodos x1 , x2 , . . . , xn en el intervalo [a, b], y
ORVFRHÀFLHQWHV c1 , c2 , . . . , cn se seleccionan para minimizar el error esperado obtenido en
la aproximación
Gauss demostró su método de b n
LQWHJUDFLyQQXPpULFDHÀFLHQWH f (x) d x ≈ ci f (xi ).
en un documento presentado a i=1
ante la Göttingen Society en
1814. Él permitió que los nodos,
DVtFRPRORVFRHÀFLHQWHVGHODV Para medir esta precisión, suponemos que la mejor elección de estos valores produce el
evaluaciones de función, fueran resultado exacto para la clase de polinomios de mayor grado, es decir, la selección que da
parámetros en la fórmula de el grado más alto de precisión.
suma y encontró la colocación /RVFRHÀFLHQWHVc1 , c2 , . . . , cn en la fórmula de aproximación son arbitrarios y los nodos
óptima de los nodos. Goldstine
[Golds], pp. 224–232, tiene una
x1 , x2 , . . . , xn están restringidos solamente por el hecho de que deben encontrarse en [a, b],
descripción interesante de este el intervalo de integración. Esto nos da 2n SDUiPHWURV D HOHJLU 6L ORV FRHÀFLHQWHV GH XQ
desarrollo. polinomio se consideran parámetros, la clase de polinomios de grado por lo menos 2n 2 1
también contiene 2n parámetros. Esto, entonces, es la mayor clase de polinomios para los
que es razonable esperar que una fórmula sea exacta. Con la elección adecuada de valores y
constantes, se puede obtener la precisión en este conjunto.
170 CAPÍTULO 4 Diferenciación numérica e integración

Para ilustrar el procedimiento para elegir los parámetros adecuados, mostraremos cómo
VHOHFFLRQDUORVFRHÀFLHQWHV\QRGRVFXDQGRn 5 2 y el intervalo de integración es [21, 1].
Entonces, analizaremos la situación más general para una selección arbitraria de nodos y
FRHÀFLHQWHV\PRVWUDUHPRVFyPRVHPRGLÀFDODWpFQLFDDOLQWHJUDUVREUHXQLQWHUYDORDUEL-
trario.
Suponga que queremos determinar c1, c2, x1, y x2 de tal forma que la fórmula de inte-
gración
1
f (x) d x ≈ c1 f (x1 ) + c2 f (x2 )
−1

da el resultado exacto siempre que f (x) es un polinomio de grado 2(2) 2 1 5 3 o menor, es

decir, cuando

f (x) = a0 + a1 x + a2 x 2 + a3 x 3 ,

para algún conjunto de constantes a0, a1, a2, y a3. Puesto que

(a0 + a1 x + a2 x 2 + a3 x 3 ) d x = a0 1 d x + a1 x d x + a2 x 2 d x + a3 x 3 d x,

esto es equivalente a mostrar que la fórmula da resultados exactos cuando f (x) es 1, x, x2, y
x3. Por lo tanto, necesitamos c1, c2, x1, y x2, de tal forma que
1 1
c1 · 1 + c 2 · 1 = 1 d x = 2, c1 · x1 + c2 · x2 = x d x = 0,
−1 −1
1 1
2
c1 · x12 + c2 · x22 = x2 dx = , y c1 · x13 + c2 · x23 = x 3 d x = 0.
−1 3 −1

Un poco de álgebra muestra que este sistema de ecuaciones tiene una solución única
√ √
3 3
c1 = 1, c2 = 1, x1 = − , y x2 = ,
3 3
lo cual da la fórmula de aproximación
√ √
1
− 3 3
f (x) d x ≈ f + f . (4.40)
−1 3 3

Esta fórmula tiene grado de precisión tres; es decir, produce el resultado exacto para cada
polinomio de grado tres o menor.

Polinomios de Legendre
/DWpFQLFDTXHKHPRVGHVFULWRVHSXHGHXWLOL]DUSDUDGHWHUPLQDUORVQRGRV\FRHÀFLHQWHVSDUD
las fórmulas que dan resultados exactos para polinomios de grado superior, pero un método
alterno los obtiene con mayor facilidad. En las secciones 8.2 y 8.3, consideraremos diferentes
conjuntos de polinomios ortogonales, funciones que tienen la propiedad de que una integral
GHÀQLGDGHOSURGXFWRGHFXDOTXLHUDGHGRVGHHOODVHV(OFRQMXQWRUHOHYDQWHSDUDQXHVWUR
problema son los polinomios de Legendre, un conjunto {P0 (x), P1 (x), . . . , Pn (x), . . . , }
con las siguientes propiedades:

(1) Para cada n, Pn (x) es un polinomio mónico de grado n.

1
(2) P(x)Pn (x) d x = 0 siempre que P(x) sea un polinomio de grado menor a n.
−1
 4.7 Cuadratura gaussiana 171

Recuerde que los polinomios Los primeros polinomios de Legendre son

mónicosWLHQHQXQFRHÀFLHQWH
principal de 1. 1
P0 (x) = 1, P1 (x) = x, P2 (x) = x 2 − ,
3
3 6 3
P3 (x) = x 3 − x, y P4 (x) = x 4 − x 2 + .
5 7 35
Adrien-Marie Legendre (1752– Las raíces de estos polinomios son distintas, se encuentran en el intervalo (21, 1), tienen
1833) introdujo este conjunto
GHSROLQRPLRVHQ7XYR
una simetría respecto al origen y, más importante, son la elección correcta para determinar
numerosas disputas prioritarias ORVSDUiPHWURVTXHQRVSURSRUFLRQDQORVQRGRV\FRHÀFLHQWHVSDUDQXHVWURPpWRGRGHFXD-
con Gauss, principalmente dratura.
debido a la falla de Gauss al Los nodos x1 , x2 , . . . , xn necesarios para producir una fórmula de aproximación integral
publicar muchos de sus resultados
originales mucho después de que
que proporcione resultados exactos para cualquier polinomio de grado menor a 2n son las
los había descubierto. raíces del polinomio de Legendre de enésimo grado. Esto se establece mediante el siguiente
resultado.

Teorema 4.7 Suponga que x1 , x2 , . . . , xn son las raíces del n-ésimo polinomio de Legendre Pn (x) y que
para cada i = 1, 2, . . . , n, los números ciHVWiQGHÀQLGRVSRU
n
1
x − xj
ci = d x.
−1 j=1 xi − x j
j =i

Si P(x) es cualquier polinomio de grado menor a 2n, entonces

1 n
P(x) d x = ci P(xi ).
−1 i=1

Demostración Consideramos primero la situación para un polinomio P(x) de grado menor a n.

Reescriba P(x HQWpUPLQRVGHORVSROLQRPLRVGHFRHÀFLHQWHVGH/DJUDQJH n 2 1)-ésimos
con nodos en las raíces del n-ésimo polinomio de Legendre Pn(x). El término de error para
esta representación implica la n-ésima derivada de P(x). Puesto que P(x) es de grado menor
a n, la n-ésima derivada de P(x) es 0 y esta representación es exacta. Por lo tanto,
n n n
x − xj
P(x) = P(xi )L i (x) = P(xi )
i=1 i=1 j=1
xi − x j
j =i

y
⎡ ⎤
n n
1 1 ⎢ x − xj ⎥
P(x) d x = ⎢ P(xi )⎥
⎣ xi − x j ⎦ dx
−1 −1 i=1 j=1
j =i
⎡ ⎤
n n n
⎢ 1
x − xj ⎥
= ⎢ dx⎥
⎣ xi − x j ⎦ P(xi ) = ci P(xi ).
i=1 −1 j=1 i=1
j =i

Por lo tanto, el resultado es verdad para polinomios de grado menor a n.

Ahora considere un polinomio P(x) de grado por lo menos n pero menor a 2n. Divida
P(x) entre el n-ésimo polinomio de Legendre Pn(x). Esto proporciona dos polinomios Q(x) y
R(x), cada uno de grado menor a n, con

P(x) = Q(x)Pn (x) + R(x).

172 CAPÍTULO 4 Diferenciación numérica e integración

Observe que xi es una raíz de Pn(x) para cada i = 1, 2, . . . , n, por lo que tenemos

P(xi ) = Q(xi )Pn (xi ) + R(xi ) = R(xi ).

Ahora, recurrimos a la potencia única de los polinomios de Legendre. En primer lugar, el

grado del polinomio Q(x) es menor a n, por lo que (mediante la propiedad (2) de Legendre),
1
Q(x)Pn (x) d x = 0.
−1

A continuación, puesto que R(x) es un polinomio de grado menor a n, el argumento de aper-

tura implica que
1 n
R(x) d x = ci R(xi ).
−1 i=1

$OMXQWDUHVWRVKHFKRVYHULÀFDPRVTXHODIyUPXODHVH[DFWDSDUDHOSROLQRPLRP(x):
1 1 1 n n
P(x) d x = [Q(x)Pn (x) + R(x)] d x = R(x) d x = ci R(xi ) = ci P(xi ).
−1 −1 −1 i=1 i=1

Las constantes ci necesarias para la regla de cuadratura se pueden generar a partir de la

ecuación en el teorema 4.7, pero tanto las constantes como las raíces de los polinomios de
Legendre se tabulan ampliamente. La tabla 4.12 enumera estos valores para n 5 2, 3, 4 y 5.

Tabla 4.12 n Raíces rn,i Coeficientes cn,i

2 0.5773502692 1.0000000000
−0.5773502692 1.0000000000
3 0.7745966692 0.5555555556
0.0000000000 0.8888888889
−0.7745966692 0.5555555556
4 0.8611363116 0.3478548451
0.3399810436 0.6521451549
−0.3399810436 0.6521451549
−0.8611363116 0.3478548451
5 0.9061798459 0.2369268850
0.5384693101 0.4786286705
0.0000000000 0.5688888889
−0.5384693101 0.4786286705
−0.9061798459 0.2369268850

1
Ejemplo 1 Aproxime e x cos x d x mediante cuadratura gaussiana con n 5 3.
−1

Solución Las entradas en la tabla 4.12 nos dan

1
e x cos x d x ≈ 0.5e0.774596692 cos 0.774596692
−1

+ 0.8 cos 0 + 0.5e−0.774596692 cos(−0.774596692)

= 1.9333904.

La integración por partes se puede utilizar para mostrar que el valor verdadero de la integral
es 1.9334214, por lo que el error absoluto es menor a 3.2 3 1025.
 4.7 Cuadratura gaussiana 173

Cuadratura gaussiana en intervalos arbitrarios

b
Una integral a f (x) d x sobre un intervalo arbitrario [a, b] se puede transformar en una
integral sobre [2@DOXWLOL]DUHOFDPELRGHYDULDEOHV YpDVHODÀJXUD 

2x − a − b 1
t= ⇐⇒ x = [(b − a)t + a + b].
b−a 2

Figura 4.17
t
(b, 1)
1
2x 2 a 2 b
t5
b2a

a b x

21
(a, 21)

Esto permite aplicar la cuadratura gaussiana a cualquier intervalo [a, b] porque

b 1
(b − a)t + (b + a) (b − a)
f (x) d x = f dt. (4.41)
a −1 2 2

3
Ejemplo 2 Considere la integral x 6 − x 2 sen(2x) d x = 317.3442466.
1
a) Compare los resultados de la fórmula cerrada de Newton-Cotes con n 5 1, la
fórmula abierta de Newton-Cotes con n 5 1, y la cuadratura gaussiana cuando n 5 2.
b) Compare los resultados de la fórmula cerrada de Newton-Cotes con n 5 2, la
fórmula abierta de Newton-Cotes con n 5 2, y la cuadratura gaussiana cuando n 5 3.

Solución a) Cada una de las fórmulas en esta parte requiere dos evaluaciones de la función
f (x) = x 6 − x 2 sen(2x). Las aproximaciones de Newton-Cotes son
2
Cerrada n = 1 : [ f (1) + f (3)] = 731.6054420;
2
3(2/3)
Abierta n = 1 : [ f (5/3) + f (7/3)] = 188.7856682.
2
La cuadratura gaussiana aplicada a este problema requiere que la integral primero se trans-
forme en un problema cuyo intervalo de integración es [21, 1]. Por medio de la ecuación
(4.41) obtenemos
3 1
x 6 − x 2 sen(2x) d x = (t + 2)6 − (t + 2)2 sen(2(t + 2)) dt.
1 −1

Entonces, la cuadratura con n 5 2 da

3
x 6 − x 2 sen(2x) d x ≈ f (−0.5773502692 + 2) + f (0.5773502692 + 2)
1

= 306.8199344.
174 CAPÍTULO 4 Diferenciación numérica e integración

b) Cada una de las fórmulas en esta parte requiere tres evaluaciones de función. Las aproxi-
maciones de Newton-Cotes son

1
Cerrada n = 2 : [ f (1) + 4 f (2) + f (3)] = 333.2380940;
3
4(1/2)
Abierta n = 2 : [2 f (1.5) − f (2) + 2 f (2.5)] = 303.5912023.
3

La cuadratura gaussiana con n 5 3, una vez que la transformación se ha realizado, obtenemos

3
x 6 − x 2 sen(2x) d x ≈ 0.5 f (−0.7745966692 + 2)
1

+ 0.8 f (2) + 0.5 f (−0.7745966692 + 2) = 317.2641516.

Los resultados de la cuadratura gaussiana son claramente superiores en cada instancia.

La sección Conjunto de ejercicios 4.7 está disponible en línea. Encuentre la ruta de

acceso en las páginas preliminares.

4.8 Integrales múltiples

/DVWpFQLFDVDQDOL]DGDVHQODVVHFFLRQHVSUHYLDVVHSXHGHQPRGLÀFDUSDUDXWLOL]DUVHHQOD
aproximación de integrales múltiples. Considere la integral doble

f (x, y) d A,
R

cuando R = { (x, y) | a ≤ x ≤ b, c ≤ y ≤ d }, para algunas constantes a, b, c y d, es

XQDUHJLyQUHFWDQJXODUHQHOSODQR 9pDVHODÀJXUD

Figura 4.18

z 5 f (x, y)

a c
d

b y
R
x

La siguiente ilustración muestra la manera en la que la regla compuesta trapezoidal por

medio de dos subintervalos en cada dirección coordinada se aplicaría a esta integral.
 4.8 Integrales múltiples 175

Ilustración Al escribir la integral doble como una integral iterada obtenemos

b d
f (x, y) d A = f (x, y) dy d x.
a c
R

3DUDVLPSOLÀFDUODQRWDFLyQVL k = (d −c)/2 y h = (b−a)/2. Aplique la regla compuesta

trapezoidal a la integral interior para obtener
d
k c+d
f (x, y) dy ≈ f (x, c) + f (x, d) + 2 f x, .
c 2 2

Esta aproximación es de orden O (d − c)3 . A continuación, aplique la regla compuesta tra-

pezoidal nuevamente para aproximar la integral de esta función de x:
b d b
d −c c+d
f (x, y) dy dx ≈ f (x, c) + 2 f x, + f (x, d) d x
a c a 4 2
b−a d −c c+d
= f (a, c) + 2 f a, + f (a, d)
4 4 2
b−a d −c a+b a+b c+d a+b
+ 2 f ,c + 2 f , + f ,d
4 4 2 2 2 2
b−a d −c c+d
+ f (b, c) + 2 f b, + f (b, d)
4 4 2
(b − a)(d − c)
= f (a, c) + f (a, d) + f (b, c) + f (b, d)
16
a+b a+b c+d c+d
+2 f ,c + f ,d + f a, + f b,
2 2 2 2
a+b c+d
+4f ,
2 2

Esta aproximación es de orden O (b − a)(d − c) (b − a)2 + (d − c)2  /D ÀJXUD

4.19 muestra una cuadrícula con el número de evaluaciones funcionales en cada uno de los
nodos utilizados en la aproximación.

Figura 4.19
y

1 2 1
d

1 2 4 2
2 (c 1 d)

1 2 1
c

a 1 x
2 (a 1 b) b
176 CAPÍTULO 4 Diferenciación numérica e integración

Como muestra la ilustración, el procedimiento es bastante sencillo. Sin embargo, el

número de evaluaciones de función crece con el cuadrado del número requerido para una
sola integral. En una situación práctica, no esperaríamos utilizar un método tan básico como
la regla compuesta trapezoidal con n 5 2. En cambio, usaremos la regla compuesta de Simp-
son que es más apropiada para ilustrar la técnica de aproximación general, a pesar de que se
podría usar cualquier otra fórmula en su lugar.
Para aplicar la regla compuesta de Simpson, dividimos la región R al subdividir [a, b]
y [c, d@HQXQQ~PHURSDUGHVXELQWHUYDORV3DUDVLPSOLÀFDUODQRWDFLyQVHOHFFLRQDPRVHQ-
teros pares n y m y subdividimos [a, b] y [c, d] con puntos de malla igualmente espaciados
x0 , x1 , . . . , xn y y0 , y1 , . . . , ym, respectivamente. Estas subdivisiones determinan tamaños
de pasos h = (b − a)/n y k = (d − c)/m. Al escribir la integral doble como la integral
iterada
b d
f (x, y) d A = f (x, y) dy d x,
a c
R

primero utilizamos la regla compuesta de Simpson para aproximar

d
f (x, y) dy,
c

tomando x como una constante.

Si y j = c + jk, para cada j = 0, 1, . . . , m. Entonces

⎡ ⎤
d (m/2)−1 m/2
k⎣
f (x, y) dy = f (x, y0 ) + 2 f (x, y2 j ) + 4 f (x, y2 j−1 ) + f (x, ym )⎦
c 3 j=1 j=1

(d − c)k 4 ∂ 4 f (x, μ)
− ,
180 ∂ y4

para algunas m en (c, d). Por lo tanto,

b d b (m/2)−1 b
k
f (x, y) d y d x = f (x, y0 ) d x + 2 f (x, y2 j ) d x
a c 3 a j=1 a

m/2 b b
+4 f (x, y2 j−1 ) d x + f (x, ym ) d x
j=1 a a

b
(d − c)k 4 ∂ 4 f (x, μ)
− d x.
180 a ∂ y4

Ahora, utilizamos la regla compuesta de Simpson para las integrales en esta ecuación. Si
xi = a + i h, para cada i = 0, 1, . . . , n. Entonces, para cada j = 0, 1, . . . , m tenemos

b (n/2)−1 n/2
h
f (x, y j ) d x = f (x0 , y j ) + 2 f (x2i , y j ) + 4 f (x2i−1 , y j ) + f (xn , y j )
a 3 i=1 i=1

(b − a)h 4 ∂ 4 f
− (ξ j , y j ),
180 ∂x4
 4.8 Integrales múltiples 177

para algunos ξ j en (a, b). La aproximación resultante tiene la forma

b d (n/2)−1
hk
f (x, y) d y d x ≈ f (x0 , y0 ) + 2 f (x2i , y0 )
a c 9 i=1
n/2
+4 f (x2i−1 , y0 ) + f (xn , y0 )
i=1
(m/2)−1 (m/2)−1 (n/2)−1
+2 f (x0 , y2 j ) + 2 f (x2i , y2 j )
j=1 j=1 i=1

(m/2)−1 n/2 (m/2)−1

+4 f (x2i−1 , y2 j ) + f (xn , y2 j )
j=1 i=1 j=1

m/2 m/2 (n/2)−1

+4 f (x0 , y2 j−1 ) + 2 f (x2i , y2 j−1 )
j=1 j=1 i=1

m/2 n/2 m/2

+4 f (x2i−1 , y2 j−1 ) + f (xn , y2 j−1 )
j=1 i=1 j=1

(n/2)−1 n/2
+ f (x0 , ym ) + 2 f (x2i , ym ) + 4 f (x2i−1 , ym )
i=1 i=1

+ f (xn , ym ) .

El término de error E está dado por

(m/2)−1 m/2
−k(b − a)h 4 ∂ 4 f (ξ0 , y0 ) ∂ 4 f (ξ2 j , y2 j ) ∂ 4 f (ξ2 j−1 , y2 j−1 )
E= +2 +4
540 ∂x4 j=1
∂x4 j=1
∂x4
b
∂ 4 f (ξm , ym ) (d − c)k 4 ∂ 4 f (x, μ)
+ − d x.
∂x4 180 a ∂ y4

Si ∂ 4 f /∂ x 4 es continua, el teorema del valor intermedio 1.11 se puede aplicar repeti-

damente para mostrar que la evaluación de las derivadas parciales respecto a x se pueden
reemplazar con un valor común y que

b
−k(b − a)h 4 ∂4 f (d − c)k 4 ∂ 4 f (x, μ)
E= 3m 4 (η, μ) − d x,
540 ∂x 180 a ∂ y4

para algunos (η, μ) en R. Si ∂ 4 f /∂ y 4 también es continua, el teorema de valor medio para

integrales implica que

b
∂ 4 f (x, μ) ∂4 f
d x = (b − a) (η̂, μ̂),
a ∂ y4 ∂ y4
178 CAPÍTULO 4 Diferenciación numérica e integración

para algunos (η̂, μ̂) en R. Puesto que m = (d − c)/k, el término de error tiene la forma

−k(b − a)h 4 ∂4 f (d − c)(b − a) 4 ∂ 4 f

E= 3m 4 (η, μ) − k (η̂, μ̂),
540 ∂x 180 ∂ y4

ORFXDOVHVLPSOLÀFDHQ

(d − c)(b − a) 4 ∂ 4 f 4∂ f
4
E =− h (η, μ) + k (η̂, μ̂) ,
180 ∂x4 ∂ y4

para algunos (η, μ) y (η̂, μ̂) en R.

Ejemplo 1 Utilice la regla compuesta de Simpson con n 5 4 y m 5 2 para aproximar

2.0 1.5
ln(x + 2y) d y d x.
1.4 1.0

Solución Los tamaños de paso para esta aplicación son h = (2.0 − 1.4)/4 = 0.15 y k =
(1.5 − 1.0)/2 = 0.25. La región de integración RVHPXHVWUDHQODÀJXUDMXQWRFRQORV
nodos (xi , y j ), donde i = 0, 1, 2, 3, 4 y j = 0, 1, 2.  7DPELpQ PXHVWUD TXH ORV FRHÀFLHQWHV
w i, j de f (xi , yi ) = ln(xi +2yi ) en la suma que da la aproximación de la regla compuesta de
Simpson para la integral.

Figura 4.20
y

1 4 2 4 1
1.50

4 16 8 16 4
1.25

1 4 2 4 1
1.00

1.40 1.55 1.70 1.85 2.00 x

La aproximación es
4 2
2.0 1.5
(0.15)(0.25)
ln(x + 2y) d y d x ≈ w i, j ln(xi + 2y j )
1.4 1.0 9 i=0 j=0

= 0.4295524387.

tenemos

∂4 f −6 ∂4 f −96
(x, y) = y (x, y) = ,
∂x 4 (x + 2y)4 ∂y 4 (x + 2y)4

y los valores máximos de los valores absolutos de estas derivadas parciales se presentan en
R cuando x 5 1.4 y y 5 1.0. Por lo que el error está acotado por

(0.5)(0.6) 6 96
|E| ≤ (0.15)4 máx + (0.25)4 máx ≤ 4.72 × 10−6 .
180 (x,y)inR (x + 2y) 4 (x,y)inR (x + 2y)4
 4.8 Integrales múltiples 179

El valor real de la integral para 10 lugares decimales es

2.0 1.5
ln(x + 2y) d y d x = 0.4295545265,
1.4 1.0

por lo que la aproximación es precisa dentro de 2.1 3 1026.

Las mismas técnicas se pueden aplicar para la aproximación de las integrales triples
así como integrales superiores para las funciones de más de tres variables. El número de
evaluaciones requeridas para la aproximación es el producto del número de evaluaciones
requeridas cuando se aplica el método a cada variable.

Cuadratura gaussiana para aproximación de integral doble

3DUDUHGXFLUHOQ~PHURGHHYDOXDFLRQHVIXQFLRQDOHVVHSXHGHQLQFOXLUPpWRGRVPiVHÀFLHQ-
tes, como la cuadratura gaussiana, la integración de Romberg o la cuadratura adaptable, en
lugar de las fórmulas de Newton-Cotes. El siguiente ejemplo ilustra el uso de cuadratura
gaussiana para la integral considerada en el ejemplo 1.

Ejemplo 2 Utilice la cuadratura gaussiana con n 5 3 en ambas dimensiones para aproximar la integral
2.0 1.5
ln(x + 2y) d y d x.
1.4 1.0

Solución Antes de emplear la cuadratura gaussiana para aproximar esta integral, necesita-
mos transformar la región de integración

R = { (x, y) | 1.4 ≤ x ≤ 2.0, 1.0 ≤ y ≤ 1.5 } en

R̂ = { (u, v) | −1 ≤ u ≤ 1, −1 ≤ v ≤ 1 }.

Las transformaciones lineales que cumplen esto son

1 1
u= (2x − 1.4 − 2.0) y v= (2y − 1.0 − 1.5),
2.0 − 1.4 1.5 − 1.0

o, de manera equivalente x = 0.3u + 1.7 y y = 0.25v + 1.25. Al emplear este cambio de

variables obtenemos una integral en la que se puede aplicar la cuadratura gaussiana:

2.0 1.5 1 1
ln(x + 2y) d y d x = 0.075 ln(0.3u + 0.5v + 4.2) dv du.
1.4 1.0 −1 −1

La fórmula de cuadratura gaussiana para n 5 3 tanto en u como en v requiere que utilicemos

los nodos

u 1 = v1 = r3,2 = 0, u 0 = v0 = r3,1 = −0.7745966692,

y

u 2 = v2 = r3,3 = 0.7745966692.

Los pesos asociados con c3,2 = 0.8 y c3,1 = c3,3 = 0.5. (Estos se muestran en la tabla 4.12
en la página 232.) La aproximación resultante es

2.0 1.5 3 3
ln(x + 2y) d y d x ≈ 0.075 c3,i c3, j ln(0.3r3,i + 0.5r3, j + 4.2)
1.4 1.0 i=1 j=1

= 0.4295545313.
180 CAPÍTULO 4 Diferenciación numérica e integración

A pesar de que este resultado sólo requiere nueve evaluaciones funcionales en comparación
con 15 para la regla compuesta de Simpson considerada en el ejemplo 1, éste es preciso den-
tro de 4.8×10−9, en comparación con la precisión 2.1 × 10−6 en el ejemplo 1.

Regiones no rectangulares
El uso de métodos de aproximación para integrales dobles no está limitado a integrales con
regiones rectangulares de integración. Las técnicas analizadas previamente se pueden modi-
ÀFDUSDUDDSUR[LPDUODVLQWHJUDOHVGREOHVGHODIRUPD
b d(x)
f (x, y) d y dx (4.42)
a c(x)

d b(y)
f (x, y) d x d y. (4.43)
c a(y)

De hecho, las integrales en las regiones que no son de este tipo, también se pueden aproximar
al realizar subdivisiones adecuadas de la región. (Consulte el ejercicio 10.)
Para describir la técnica relacionada con la aproximación de una integral de la forma
b d(x)
f (x, y) d y d x,
a c(x)

utilizaremos la regla básica de Simpson para integrar respecto a ambas variables. El tamaño
de paso para la variable x es h = (b − a)/2, pero el tamaño de paso para y varía con x (véase
ODÀJXUD \VHHVFULEH

d(x) − c(x)
k(x) = .
2

Figura 4.21

z 5 f (x, y)

d(a) y 5 d(x) A(x)

d(b)
k(a)
k(b)
c(b) a
c(a)
y 5 c(x)
k(a 1 h) y
b
R y 5 d(x)
a a1h b x
x y 5 c(x)
a) b)
 4.8 Integrales múltiples 181

Esto nos da
b d(x) b
k(x)
f (x, y) d y d x ≈ [ f (x, c(x)) + 4 f (x, c(x) + k(x)) + f (x, d(x))] d x
a c(x) a 3
h k(a)
≈ [ f (a, c(a)) + 4 f (a, c(a) + k(a)) + f (a, d(a))]
3 3
4k(a + h)
+ [ f (a + h, c(a + h)) + 4 f (a + h, c(a + h)
3
+ k(a + h)) + f (a + h, d(a + h))]
k(b)
+ [ f (b, c(b)) + 4 f (b, c(b) + k(b)) + f (b, d(b))] .
3

El algoritmo 4.4. aplica la regla compuesta de Simpson para una integral de la forma (4.42).
Las integrales de la forma (4.43) pueden, por supuesto, manejarse de manera similar.

ALGORITMO Integral doble de Simpson

4.4 Para aproximar la integral
b d(x)
I = f (x, y) d y d x :
a c(x)

ENTRADA extremos a, b: enteros positivos pares m, n.

SALIDA aproximación J para I.
Paso 1 Haga h = (b − a)/n;
J1 = 0; (Términos finales.)
J2 = 0; (Términos pares.)
J3 = 0. (Términos impares.)
Paso 2 Para i = 0, 1, . . . , n haga los pasos 3–8.
Paso 3 Haga x = a + i h; (Método compuesto de Simpson para x.)
HX = (d(x) − c(x))/m;
K 1 = f (x, c(x)) + f (x, d(x)); (Términos finales.)
K 2 = 0; (Términos pares.)
K 3 = 0. (Términos impares.)
Paso 4 Para j = 1, 2, . . . , m − 1 haga los pasos 5 y 6.
Paso 5 Haga y = c(x) + jHX;
Q = f (x, y).
Paso 6 Si j incluso determine entonces K 2 = K 2 + Q
también determine K 3 = K 3 + Q.
Paso 7 Haga L = (K 1 + 2K 2 + 4K 3 )HX/3.
d(xi )
L≈ f (xi , y) dy mediante el método compuesto de Simpson .
c(xi )

Paso 8 Si i = 0 o i = n entonces detemine J1 = J1 + L

también si i incluso determine entonces J2 = J2 + L
también determine J3 = J3 + L. (Paso final 2)
Paso 9 Haga J = h(J1 + 2J2 + 4J3 )/3.
Paso 10 SALIDA ( J );
PARE.
182 CAPÍTULO 4 Diferenciación numérica e integración

Aplicando la cuadratura gaussiana a la integral doble

b d(x)
f (x, y) dy dx
a c(x)

primero debemos transformar, para cada x en [a, b], la variable y en el intervalo [c(x), d(x)]
en la variable t en el intervalo [21, 1]. Esta transformación lineal nos da

(d(x) − c(x))t + d(x) + c(x) d(x) − c(x)

f (x, y) = f x, y dy = dt.
2 2

El cálculo reducido hace que, en Entonces, para cada x en [a, b], aplicamos la cuadratura gaussiana a la integral resultante
general, valga la pena aplicar la
cuadratura gaussiana en lugar d(x)
de una técnica de Simpson al
1
(d(x) − c(x))t + d(x) + c(x)
f (x, y) dy = f x, dt
aproximar las integrales dobles. c(x) −1 2

para producir
b d(x)
f (x, y) d y d x
a c(x)
b n
d(x) − c(x) (d(x) − c(x))rn, j + d(x) + c(x)
≈ cn, j f x, d x,
a 2 j=1
2

mientras, como antes, las raíces rn,j\ORVFRHÀFLHQWHVcn,j provienen de la tabla 4.12 en la

página 172. Ahora, el intervalo [a, b] se transforma en [21, 1], y la cuadratura gaussiana se
aplica para aproximar la integral en el lado derecho de esta ecuación. Los detalles se inclu-
yen en el algoritmo 4.5.

ALGORITMO Integral doble gaussiana

4.5 Para aproximar la integral
b d(x)
f (x, y) d y d x :
a c(x)

ENTRADA extremos a, b: enteros positivos m, n.

(Las raíces ri, j y los coeficientes c i, j necesitan estar disponibles para i = máx{m, n}
y para 1 ≤ j ≤ i.)
SALIDA aproximación J para I.
Paso 1 Haga h 1 = (b − a)/2;
h 2 = (b + a)/2;
J = 0.
Paso 2 Para i = 1, 2, . . . , m haga los pasos 3–5.
Paso 3 Haga J X = 0;
x = h 1rm,i + h 2 ;
d1 = d(x);
c1 = c(x);
k1 = (d1 − c1 )/2;
k2 = (d1 + c1 )/2.
Paso 4 Para j = 1, 2, . . . , n haga
y = k1rn, j + k2 ;
Q = f (x, y);
J X = J X + cn, j Q.
 4.8 Integrales múltiples 183

Paso 5 Haga J = J + cm,i k1 J X . (Paso final 2.)

Paso 6 Haga J = h 1 J .
Paso 7 SALIDA (J );
PARE.

Ilustración (OYROXPHQGHOVyOLGRHQODÀJXUDVHDSUR[LPDDSOLFDQGRHODOJRULWPRGHODLQWHJUDO
doble de Simpson con n 5 m 5 10 a
0.5 x2
e y/x d y d x.
0.1 x3

Esto requiere 121 evaluaciones de la función f (x, y) = e y/x y produce el valor

HOFXDODSUR[LPDHOYROXPHQGHOVyOLGRPRVWUDGRHQODÀJXUDDDSUR[LPDGD-
mente siete lugares decimales. Al aplicar el algoritmo de cuadratura gaussiana con n 5 m 5 5
sólo se requieren 25 evaluaciones de función y da la aproximación 0.03330556611, lo cual
es preciso para 11 lugares decimales.

Figura 4.22
z

1 (0.1, 0.01, e0.1) (0.5, 0.25, e0.5)

(0.1, 0.001, e0.01)

y

0.25 (0.5, 0.125, e0.25)

0.125

0.1

R (0.5, 0.25, 0)

(0.5, 0.125, 0)
0.5
x

Aproximación de integral triple

Las integrales triples de la forma

El cálculo reducido casi siempre b d(x) β(x,y)

hace que valga la pena aplicar f (x, y, z) dz dy dx
la cuadratura gaussiana en lugar a c(x) α(x,y)
de la técnica de Simpson al
aproximar integrales triples o
superiores. 9pDVHODÀJXUD VHDSUR[LPDQGHPDQHUDVLPLODU'HELGRDOQ~PHURGHFiOFXORVLPSOL-
cados, la cuadratura gaussiana es el método de elección. El algoritmo 4.6 implementa este
procedimiento.
184 CAPÍTULO 4 Diferenciación numérica e integración

Figura 4.23
z

z 5 b(x, y)

z 5 a(x, y)

a
x
y
b y 5 c(x)
R y 5 d(x)
x

ALGORITMO Integral triple gaussiana

4.6 Para aproximar la integral
b d(x) β(x,y)
f (x, y, z) dz dy d x :
a c(x) α(x,y)

ENTRADA extremos a, b; enteros positivos m, n, p.

(Las raíces ri, j y los coeficientes ci, j necesitan estar disponibles para i = máx{n, m, p}
y para 1 ≤ j ≤ i.)
SALIDA aproximación J para I.
Paso 1 Haga h 1 = (b − a)/2;
h 2 = (b + a)/2;
J = 0.
Paso 2 Para i = 1, 2, . . . , m haga los pasos 3–8.
Paso 3 Haga JX = 0;
x = h 1rm,i + h 2 ;
d1 = d(x);
c1 = c(x);
k1 = (d1 − c1 )/2;
k2 = (d1 + c1 )/2.
Paso 4 Para j = 1, 2, . . . , n haga los pasos 5–7.
Paso 5 Haga JY = 0;
y = k1rn, j + k2 ;
β1 = β(x, y);
α1 = α(x, y);
l1 = (β1 − α1 )/2;
l2 = (β1 + α1 )/2.
 4.8 Integrales múltiples 185

Paso 6 Para k = 1, 2, . . . , p haga

Haga z = l1r p,k + l2 ;
Q = f (x, y, z);
JY = JY + c p,k Q.
Paso 7 Haga JX = JX + cn, j l1 JY. (Fin del paso 4)
Paso 8 Haga J = J + cm,i k1 JX. (Fin del paso 2)
Paso 9 Haga J = h 1 J .
Paso 10 SALIDA (J );
PARE.

El siguiente ejemplo requiere la evaluación de cuatro integrales triples.

Ilustración El centro de masa de una región sólida DFRQIXQFLyQGHGHQVLGDGƱVHSUHVHQWDHQ

M yz Mx z Mx y
(x, y, z) = , , ,
M M M

donde

M yz = xσ (x, y, z) d V, Mx z = yσ (x, y, z) d V
D D

Mx y = zσ (x, y, z) d V
D

son los momentos alrededor de los planos coordenados y la masa de D es

M= σ (x, y, z) d V.
D

(OVyOLGRPRVWUDGRHQODÀJXUDHVWiDFRWDGRSRUODSDUWHVXSHULRUGHOFRQRz 2 = x 2 + y 2
y el plano z 5 2. Suponga que este sólido tiene una función de densidad dada por

σ (x, y, z) = x 2 + y2.

Figura 4.24
z

1
2 1
x 2
y
186 CAPÍTULO 4 Diferenciación numérica e integración

Al aplicar el algoritmo de integral triple gaussiana 4.6 con n 5 m 5 p 5 5 requiere 125 eva-
luaciones de función por integral y se obtienen las siguientes aproximaciones
√
2 4−x 2 2
M= √ √ x 2 + y 2 dz dy d x
−2 − 4−x 2 x 2 +y 2
√
2 4−x 2 2
=4 √ x 2 + y 2 dz dy d x ≈ 8.37504476,
0 0 x 2 +y 2
√
2 4−x 2 2
M yz = √ √ x x 2 + y 2 dz dy d x ≈ −5.55111512 × 10−17 ,
−2 − 4−x 2 x 2 +y 2
√
2 4−x 2 2
Mx z = √ √ y x 2 + y 2 dz dy d x ≈ −8.01513675 × 10−17 y
−2 − 4−x 2 x 2 +y 2
√
2 4−x 2 2
Mx y = √ √ z x 2 + y 2 dz dy d x ≈ 13.40038156.
−2 − 4−x 2 x 2 +y 2

Esto implica que la ubicación aproximada del centro de la masa es

(x, y, z) = (0, 0, 1.60003701).

Estas integrales son bastante fáciles de evaluar de manera directa. Si lo hace, descubrirá que
el centro de masa exacto se presenta en (0, 0, 1.6).

La sección Conjunto de ejercicios 4.8 está disponible en línea. Encuentre la ruta de

acceso en las páginas preliminares.

4.9 Integrales impropias

Las integrales impropias resultan cuando la noción de integración se amplía, ya sea en un
intervalo de integración en el que la función no está acotada en un intervalo con uno o más
H[WUHPRVLQÀQLWRV(QFXDOTXLHUFLUFXQVWDQFLDODVUHJODVQRUPDOHVGHXQDDSUR[LPDFLyQLQ-
WHJUDOVHGHEHQPRGLÀFDU

Singularidad del extremo izquierdo

Primero consideraremos la situación en la que el integrando no está limitado en el extremo
L]TXLHUGRGHOLQWHUYDORGHLQWHJUDFLyQFRPRVHPXHVWUDHQODÀJXUD(QHVWHFDVRGH-
cimos que f tiene una singularidad en el extremo a. A continuación, mostramos cómo se
pueden reducir las integrales impropias para problemas de este tipo.
Figura 4.25
y

y 5 f(x)

a b x
 4.9 Integrales impropias 187

En cálculo se muestra que la integral impropia con una singularidad en el extremo iz-
quierdo,
b
dx
p
,
a (x − a)

converge si y sólo si 0 < p < 1, \HQHVWHFDVRGHÀQLPRV

b x=b
1 (x − a)1− p (b − a)1− p
d x = lím = .
a (x − a) p M→a + 1− p x=M 1− p

1 1
1 1
Ejemplo 1 Muestre que la integral impropia √ d x converge pero que d x diverge.
0 x 0 x2
Solución Para la primera integral, tenemos

1 1
1 x=1
√ d x = lím x −1/2 d x = lím 2x 1/2 x=M
= 2 − 0 = 2,
0 x M→0+ M M→0+

pero la segunda integral

1 1
1 x=1
d x = lím x −2 d x = lím −x −1 x=M
0 x2 M→0+ M M→0+

HVLQÀQLWD

Si f es una función que se puede escribir en la forma

g(x)
f (x) = ,
(x − a) p

donde 0 < p < 1 y g es continua en [a, b], entonces la integral impropia

b
f (x) d x
a

también existe. Nosotros aproximaremos esta integral mediante la regla compuesta de Simp-
son, siempre que g ∈ C 5 [a, b].  (Q HVH FDVR FRQVWUXLPRV HO FXDUWR SROLQRPLR GH 7D\ORU
P4(x), para g alrededor de a,

g (a) g (a) g (4) (a)

P4 (x) = g(a) + g (a)(x − a) + (x − a)2 + (x − a)3 + (x − a)4 ,
2! 3! 4!

y escribimos
b b b
g(x) − P4 (x) P4 (x)
f (x) d x = dx + d x. (4.44)
a a (x − a) p a (x − a) p

Puesto que P(x) es un polinomio, podemos determinar con exactitud el valor de

b b
g (k) (a) g (k) (a)
4 4
P4 (x)
dx = (x − a)k− p d x = (b − a)k+1− p .
a (x − a) p k=0 a k! k=0
k!(k + 1 − p)
(4.45)
188 CAPÍTULO 4 Diferenciación numérica e integración

Por lo general, la parte dominante de la aproximación, en especial cuando el polinomio de

7D\ORUP4(x) concuerda de cerca con g(x) en todo el intervalo [a, b].
Para aproximar la integral de f, debemos añadir este valor a la aproximación de
b
g(x) − P4 (x)
d x.
a (x − a) p

3DUDGHWHUPLQDUHVWRSULPHURGHÀQLPRV
g(x)−P4 (x)
(x−a) p
, si a < x ≤ b,
G(x) =
0, si x = a.

(k)
Esto nos da una función continua en [a, b]. De hecho, 0 < p < 1 y P4 (a) concuerda con
g (k) (a) para cada k = 0, 1, 2, 3, 4, por lo que tenemos G ∈ C 4 [a, b]. Esto implica que la
regla compuesta de Simpson se puede aplicar para aproximar la integral de G sobre [a, b].
Al añadir esta aproximación al valor de la ecuación (4.45) obtenemos una aproximación
para la integral impropia de f en [a, b], dentro de la precisión de la aproximación de la regla
compuesta de Simpson.

Ejemplo 2 Utilice la regla compuesta de Simpson con h 5 0.25 para aproximar el valor de la integral
impropia
1
ex
√ d x.
0 x

Solución (OFXDUWRSROLQRPLRGH7D\ORUSDUDex alrededor de x 5 0 es

x2 x3 x4
P4 (x) = 1 + x + + + ,
2 6 24
1
ex
por lo que la parte dominante de la aproximación para √ d x es
0 x
1
P4 (x) 1
1 1 1
√ dx = x −1/2 + x 1/2 + x 3/2 + x 5/2 + x 7/2 dx
0 x 0 2 6 24
1
2 1 1 1 9/2
= lím 2x 1/2 + x 3/2 + x 5/2 + x 7/2 + x
M→0+ 3 5 21 108 M
2 1 1 1
=2+ + + + ≈ 2.9235450.
3 5 21 108
1
ex
Para la segunda parte de la aproximación para √ d x , necesitamos aproximar
0 x
1
G(x) d x , donde
0
⎧
⎨ √1 e x − P (x) , si 0 < x ≤ 1,
4
G(x) = x
⎩
0, si x = 0.
Tabla 4.13
x G(x) La tabla 4.13 enumera los valores necesarios para la regla compuesta de Simpson para esta
aproximación.
0.00 0 Por medio de estos datos y la regla compuesta de Simpson obtenemos
0.25 0.0000170
1
0.50 0.0004013 0.25
G(x) d x ≈ [0 + 4(0.0000170) + 2(0.0004013) + 4(0.0026026) + 0.0099485]
0.75 0.0026026 0 3
1.00 0.0099485
= 0.0017691.
 4.9 Integrales impropias 189

Por lo tanto,
1
ex
√ d x ≈ 2.9235450 + 0.0017691 = 2.9253141.
0 x

Este resultado es preciso dentro de la precisión de la aproximación de la regla compuesta

de Simpson para la función G. Puesto que |G (4) (x)| < 1 en [0, 1], el error está acotado por

1−0
(0.25)4 = 0.0000217.
180

Singularidad en el extremo derecho

Para aproximar la integral impropia con una singularidad en el extremo derecho, podemos
desarrollar una técnica similar, pero expandiendo los términos del extremo derecho b en lu-
gar del extremo izquierdo a. De forma alternativa, podemos realizar la sustitución

z = −x, dz = − d x
para cambiar la integral impropia en una de la forma
b −a
f (x) d x = f (−z) dz, (4.46)
a −b

la cual tiene su singularidad en el extremo izquierdo. Entonces, podemos aplicar la técnica

GHVLQJXODULGDGGHOH[WUHPRL]TXLHUGRTXH\DKDEtDPRVGHVDUUROODGR 9pDVHODÀJXUD

Figura 4.26
y Para z 5 2x y

y 5 f (x) y 5 f (2z)

a b x 2b 2a z

Una integral impropia con una singularidad en c, donde a < c < b, se trata como la
suma de las integrales impropias con singularidades en los extremos ya que
b c b
f (x) d x = f (x) d x + f (x) d x.
a a c

Singularidad infinita
(ORWURWLSRGHLQWHJUDOLPSURSLDLPSOLFDOtPLWHVLQÀQLWRVGHLQWHJUDFLyQ/DLQWHJUDOEiVLFD
de este tipo tiene la forma
∞
1
d x,
a xp
190 CAPÍTULO 4 Diferenciación numérica e integración

para p > 1. Esto se convierte en una integral con singularidad de extremo izquierdo en 0 al
realizar la sustitución de integración

t = x −1 , dt = −x −2 d x, por lo que d x = −x 2 dt = −t −2 dt.

Entonces
∞
1 0
tp 1/a
1
dx = − dt = dt.
a xp 1/a t2 0 t 2− p

De forma similar, el cambio de variable t = x −1 convierte la integral impropia

∞
a f (x) d x en una que tiene singularidad de extremo izquierdo en cero:
∞ 1/a
1
f (x) d x = t −2 f dt. (4.47)
a 0 t

Ahora, se puede aproximar por medio de la fórmula de cuadratura del tipo descrito anterior-
mente.

Ejemplo 3 Aproxime el valor de la integral impropia

∞
1
I = x −3/2 sen d x.
1 x

Solución Primero realizamos el cambio de variable t = x −1, el cual convierte la singulari-

GDGLQÀQLWDHQXQDFRQXQDVLQJXODULGDGGHH[WUHPRL]TXLHUGR(QWRQFHV
1
dt = −x −2 d x, por lo que d x = −x 2 dt = − dt,
t2
y
x=∞ t=0 −3/2 1
1 1 1
I = x −3/2 sen dx = sen t − dt = t −1/2 sen t dt.
x=1 x t=1 t t2 0

(OFXDUWRSROLQRPLRGH7D\ORUP4(t), para sen t alrededor de 0 es

1
P4 (t) = t − t 3 ,
6
por lo tanto, ⎧
⎪
⎨ sen t − t + 6 t ,
1 3
si 0 < t ≤ 1
G(t) = t 1/2
⎪
⎩0, si t = 0

está en C 4 [0, 1], y tenemos

1
1 1 sen t − t + 16 t 3
I = t −1/2 t − t 3 dt + dt
0 6 0 t 1/2
2 3/2 1 1 1 sen t − t + 16 t 3
= t − t 7/2 + dt
3 21 0 0 t 1/2
1 sen t − t + 16 t 3
= 0.61904761 + dt.
0 t 1/2
 4.10 Software numérico y revisión del capítulo 191

El resultado a partir de la regla compuesta de Simpson con n 5 16 para la integral restante

HV(VWRGDXQDDSUR[LPDFLyQÀQDOGH

I = 0.0014890097 + 0.61904761 = 0.62053661,

lo cual es preciso dentro de 4.0 × 10−8 .

La sección Conjunto de ejercicios 4.9 está disponible en línea. Encuentre la ruta de

acceso en las páginas preliminares.

4.10 Software numérico y revisión del capítulo

La mayoría del software para integrar una función de una única variable real está basado
ya sea en el enfoque adaptable o en fórmulas gaussianas extremadamente precisas. La in-
WHJUDFLyQFDXWHORVDGH5RPEHUJHVXQDWpFQLFDDGDSWDEOHTXHLQFOX\HXQDYHULÀFDFLyQSDUD
garantizar que el integrando se comporta suavemente sobre los subintervalos de la integral
de integración. Este método se ha utilizado con éxito en bibliotecas de software. En general,
múltiples integrales se aproximan con la ampliación de buenos métodos adaptables hasta
dimensiones superiores. La cuadratura del tipo gaussiano también se recomienda para dis-
minuir el número de evaluaciones de función.
Las rutinas principales en las bibliotecas IMSL y NAG están basadas en QUADPACK:
A Subroutine Package for Automatic Integration (QUADPACK: Un paquete de subrutinas
para integración automática) de R. Piessens, E. de Doncker-Kapenga, C. W. Uberhuber y D.
K. Kahaner y publicado por Springer-Verlag en 1983 [PDUK].
La biblioteca IMSL contiene un esquema de integración adaptable con base en la regla
Gaussiana-Kronrod de 21 puntos mediante la regla gaussiana de 10 puntos para cálculo
de error. Las reglas gaussianas utilizan 10 puntos x1 , . . . , x10 y pesos w 1 , . . . , w 10 para la
10 b
fórmula de cuadratura i=1 w i f (xi ) para aproximar a f (x) d x . A continuación se utilizan
los puntos adicionales x11 , . . . , x21, y los pesos nuevos v1 , . . . , v21, en la fórmula de Kron-
21
rod i=1 vi f (xi ). Los resultados de las dos fórmulas se comparan para eliminar el error.
La ventaja de utilizar x1 , . . . , x10 en cada fórmula es que f sólo necesita evaluarse en 21
puntos. Si se utilizaran reglas gaussianas independientes de 10 y 21 puntos, se necesitarían
31 evaluaciones de función. Este procedimiento permite singularidades de extremo en el
integrando.
2WUDVVXEUXWLQDV,06/SHUPLWHQVLQJXODULGDGHVGHH[WUHPRVLQJXODULGDGHVHVSHFLÀFDGDV
SRUHOXVXDULRHLQWHUYDORVLQÀQLWRVGHLQWHJUDFLyQ$GHPiVH[LVWHQUXWLQDVSDUDDSOLFDUUHJODVGH
Gauss-Kronrod para integrar una función de dos variables y una rutina para utilizar cuadratura
para integrar una función de n variables sobre n intervalos de la forma [ai, bi].
La Biblioteca NAG incluye una rutina para calcular la integral de f sobre el intervalo
[a, b] mediante un método adaptable con base en cuadratura gaussiana mediante reglas de
.URQURGGHSXQWRV\GH*DXVVGHSXQWRV7DPELpQWLHQHXQDUXWLQDSDUDDSUR[LPDU
una integral mediante una familia de fórmulas tipo gaussianas con base en 1, 3, 5, 7, 15, 31,
63, 127 y 255 nodos. Estas reglas entrelazadas de alta precisión se deben a Patterson [Pat]
y se utilizan de manera adaptable. NAG incluye muchas otras subrutinas para aproximar
integrales.
A pesar de que la diferenciación numérica es inestable, se necesitan fórmulas de aproxi-
mación de derivadas para resolver ecuaciones diferenciales. La Biblioteca NAG incluye una
subrutina para la diferenciación numérica de una función de una variable real con diferencia-
ción para que la catorceava derivada sea posible. IMSL tiene una función que usa un cambio
DGDSWDEOHHQWDPDxRGHSDVRSDUDGLIHUHQFLDVÀQLWDVSDUDDSUR[LPDUODSULPHUDVHJXQGDR
tercera derivada de f en x dentro de una tolerancia determinada. IMSL también incluye una
VXEUXWLQDSDUDFDOFXODUODVGHULYDGDVGHXQDIXQFLyQGHÀQLGDHQXQFRQMXQWRGHSXQWRVPH-
192 CAPÍTULO 4 Diferenciación numérica e integración

diante interpolación cuadrática. Ambos paquetes permiten la diferenciación e integración de

splines cúbicos interpolantes construidos por las subrutinas mencionadas en la sección 3.5.

Las secciones Preguntas de análisis, Conceptos clave y Revisión del capítulo están dis-
ponibles en línea. Encuentre la ruta de acceso en las páginas preliminares.