trapecios y simpson

Vista previa del material en texto

88 CÁLCULO / CIENCIAS AMBIENTALES / TEMA 5
5.5 Integración numérica
Métodos de Newton-Côtes
De cara a calcular la integral definida:∫ b
a
f(x) dx
se llaman Métodos de Newton-Côtes a los que se basan en integrar, en lugar de la función
dada f(x), un polinomio de interpolación que aproxime a f(x) en [a, b]. Se trata por
tanto de toda una familia general de métodos, según el polinomio de interpolación que se
considere (puede elegirse diferente grado, diferentes puntos para interpolar, etc.). Para
el caso de las interpolaciones lineal y cuadrática, estos métodos se denominan Método
de los Trapecios y Método de Simpson, respectivamente.
Método de los trapecios
Como se ha comentado, el Método de los trapecios es un Método de Newton-Côtes
basado en la interpolación lineal.
La idea esencial por tanto, de cara a integrar f(x) desde el punto (a, f(a)) hasta
(b, f(b)), es aproximar f(x) por su polinomio de interpolación lineal en [a, b] (ver figura).
f(x) ≈ P1(x) =
x− b
a− b
f(a) +
x− a
b− a
f(b) , ∀x ∈ [a, b]
y aśı:
I =
∫ b
a
f(x) dx ≃
∫ b
a
P1(x) dx =
b− a
2
(f(a) + f(b))
a b
x
fHxL
a b
x
P1HxL
En definitiva se trata de aproximar el valor de la integral I por el área del trapecio
que determinan las rectas x = a, x = b, el eje de abscisas y la recta que une los puntos:
(a, f(a)) y (b, f(b)).
CÁLCULO / CIENCIAS AMBIENTALES / TEMA 5 89
Si recordamos la expresión del error de la interpolación lineal, suponiendo que f(x)
es continua y derivable dos veces en el intervalo [a, b]:
f(x) = P1(x) + ε(x)
ε(x) =
f ′′(ξ)
2
(x− a)(x− b), a ≤ ξ ≤ b
Tendremos entonces que:
I =
∫ b
a
f(x)dx =
b− a
2
(f(a) + f(b)) + E
donde el error de la integración numérica E será, obviamente:
E =
∫ b
a
ε(x)dx =
f ′′(ξ)
2
∫ b
a
(x− a)(x− b) dx
Integrando en esta última expresión y denominando h = b−a se concluye fácilmente
en que:
E = −h
3
12
f ′′(ξ) ⇒ |E| ≤
∣∣∣∣h312 M2
∣∣∣∣
siendo M2 el valor máximo que alcance la derivada segunda de la función en el intervalo
dado [a, b].
Método de los Trapecios compuesto
Si el intervalo en el que se realiza la integral es grande, el Método de los Trapecios Simple
suele ser muy impreciso. Para mejorar la exactitud, es posible subdividir el intervalo en
otros más pequeños y aplicar en cada uno de ellos el Método simple.
De esta manera, el Método de los Trapecios compuesto o generalizado consiste en
tomar una partición P = {x0, x1, . . . , xn} de [a, b], (x0 = a, xn = b), equiespaciada, es
decir: xi+1 − xi = h, ∀i = 1, . . . , n. Tendremos aśı que:
h =
b− a
n
Teniendo en cuenta las propiedades básicas de la integral definida:∫ b
a
f(x) dx =
∫ x1
x0
f(x)dx+
∫ x2
x1
f(x)dx+ . . .+
∫ xn
xn−1
f(x)dx
y aplicando a cada integral el Método simple:∫ b
a
f(x) dx ≈ h
2
(f(x0) + f(x1)) +
h
2
(f(x1) + f(x2)) + . . .+
h
2
(f(xn−1) + f(xn)) =
=
h
2
(f(x0) + 2 (f(x1) + f(x2) + . . .+ f(xn−1)) + f(xn))
90 CÁLCULO / CIENCIAS AMBIENTALES / TEMA 5
Tenemos por tanto la expresión final para el Método de los Trapecios Generalizado:∫ b
a
f(x) dx ≈ h
2
(
f(a) + 2
n−1∑
i=1
f(xi) + f(b)
)
En lo que respecta al error de integración, será evidentemente igual a la suma de los
errores de cada una de las aplicaciones del método simple:
E = E1 + E2 + . . .+ En = −
h3
12
f ′′(ξ1)−
h3
12
f ′′(ξ2)− . . .−
h2
12
f ′′(ξn)
si denominamos M2 al máximo de la función f
′′(x) en [a, b] tendremos finalmente:
|E| ≤
∣∣∣∣h312nM2
∣∣∣∣ = ∣∣∣∣(b− a)12 h2M2
∣∣∣∣
Tomaremos habitualmente E definido no negativo, por lo que es frecuente escribir
directamente:
E ≤
∣∣∣∣h312nM2
∣∣∣∣ = ∣∣∣∣(b− a)12 h2M2
∣∣∣∣
obviando el valor absoluto para E.
Ejemplo: Calcular el valor aproximado de la integral,∫ 1
0
xdx
(x+ 1)(x+ 2)
utilizando la regla de los trapecios compuesta con n = 8 subintervalos. Evaluar exactamente el
valor de la integral y compárese con el valor aproximado obtenido.
De forma exacta:
I =
∫ 1
0
x
(x+ 1)(x+ 2)
dx =
x
(x+ 1)(x+ 2)
=
A
x+ 1
+
B
x+ 2
=
A(x+ 2) +B(x+ 1)
(x+ 1)(x+ 2)
⇒ x = A(x+ 2) +B(x+ 1) ⇒
{
x = −1 ⇒ A = −1
x = −2 ⇒ B = 2
=
∫ 1
0
(
−1
x+ 1
+
2
x+ 2
)
dx = − log(x+ 1) + 2 log(x+ 2)|10 = log
(x+ 2)2
(x+ 1)
∣∣∣∣1
0
=
= log
9
2
− log 4 = 0.1177830
Método de los Trapecios, con n = 8.
Dividimos el intervalo [0, 1] en 8 subintervalos y calculamos los correspondientes valores del
integrando:
CÁLCULO / CIENCIAS AMBIENTALES / TEMA 5 91
x0 x1 x2 x3 x4 x5 x6 x7 x8
0. 0.125 0.25 0.375 0.5 0.625 0.75 0.875 1.0
f(x0) f(x1) f(x2) f(x3) f(x4) f(x5) f(x6) f(x7) f(x8)
0. 0.05228 0.08888 0.11483 0.13333 0.14652 0.15584 0.162319 0.16666
Finalmente, aplicamos la fórmula antes deducida:
I ≈ h
2
[f(x0) + 2 (f(x1) + f(x2) + f(x3) + f(x4) + f(x5) + f(x6) + f(x7)) + f(x8)]
≈ 0.125
2
[0 + 2(0.05228 + 0.0888 + 0.11483 + 0.1333 + 0.14652 + 0.15584 + 0.162319) + 0.1666]
≈ 0.117166
que da una buena aproximación al resultado exacto. En la próxima sección completaremos este
ejercicio mediante el uso del Método de Simpson y comprobaremos que proporciona una mejor
aún aproximación.
Si realizamos el mismo cálculo con un número diferente de subintervalos, se obtienen los
siguientes resultados:
n In
n = 1 0.08333
n = 2 0.108333
n = 3 0.113492
n = 4 0.11535
n = 5 0.11622
n = 10 0.11739
n = 100 0.1177791
Método de Simpson
El Método de Simpson es un método de Newton-Côtes de segundo orden, es decir basado
en integrar un polinomio de interpolación de segundo grado, de la forma siguiente:
Dada la función f(x) en [a, b], tomaremos como tercer punto para la interpolación
el punto medio de dicho intervalo, es decir: xm =
a+b
2 , y denominaremos h =
b−a
2 a la
semianchura del intervalo. De esta forma el polinomio de interpolación de grado 2 que
pasa por (a, f(a)), (xm, f(xm)) y (b, f(b)) será:
P2(x) = f(a) +
f(xm)− f(a)
h
(x− a) + f(a) + f(b)− 2f(xm)
2h2
(x− a)(x− xm)
No es dif́ıcil calcular la integral de P2(x) entre a y b, de manera que se obtiene:∫ b
a
f(x) dx ≈
∫ b
a
P2(x) dx =
h
3
(f(a) + 4f(xm) + f(b))
92 CÁLCULO / CIENCIAS AMBIENTALES / TEMA 5
fórmula del Método de Simpson (o Método de Simpson simple).
La evaluación del error de integración da lugar a un curioso resultado.
Suponiendo que la función f(x) es derivable al menos cuatro veces en el intervalo considerado,
podemos desarrollar por la fórmula de Taylor la función f(x) en x = xm hasta tercer orden (resto
de Taylor de orden 4):
f(x) = P3(x) +R4(x) =
= f(xm) + f
′(xm)(x− xm) +
f ′′(xm)
2
(x− xm)2 +
f ′′′(xm)
3!
(x− xm)3 +R4(x)
con
R4(x) =
f (4)(ξ)
4!
(x− xm)4
De esta manera tendremos:
f(a) = f(xm − h) = f(xm) + f ′(xm)(−h) +
f ′′(xm)
2
(−h)2 + f
′′′(xm)
3!
(−h)3 + f
(4)(ξ)
4!
(−h)4
f(b) = f(xm + h) = f(xm) + f
′(xm)h+
f ′′(xm)
2
h2 +
f ′′′(xm)
3!
h3 +
f (4)(ξ)
4!
h4
Con un breve cálculo se concluye en la expresión (para la fórmula del Método de Simpson):
h
3
(f(a) + 4f(xm) + f(b)) =
h
3
(
6f(xm) + f
′′(xm)h
2 +
1
12
f (4)(ξ)h4
)
=
= 2hf(xm) +
f ′′(xm)
3
h3 +
1
36
f (4)(ξ)h5
Por otro lado, si integramos el desarrollo de Taylor tendremos (simplificando los resultados):∫ b
a
f(x)dx =
∫ b
a
(P3(x) +R4(x)) dx =
=
∫ b
a
(
f(xm) + f
′(xm)(x− xm) +
f ′′(xm)
2
(x− xm)2 +
f ′′′(xm)
3!
(x− xm)3 +R4(x)
)
dx =
= 2hf(xm) +
f ′′(xm)
3
h3 +
f (4)(ξ)
60
h5
Finalmente el error de integración no es más que (tomando nuevamente el error como definido
positivo):
E =
∣∣∣∣∣
∫ b
a
f(x)dx − h
3
(f(a) + 4f(xm) + f(b))
∣∣∣∣∣
de manera que:
E =
∣∣∣∣f (4)(ξ)60 h5 − f (4)(ξ)36 h5
∣∣∣∣ = ∣∣∣∣− 190 f (4)(ξ)h5
∣∣∣∣
Si denominamos M4 al máximo que alcance la derivada cuarta de la función en el intervalo
[a, b], tendremos finalmente:
E ≤
∣∣∣∣ 190 h5 M4
∣∣∣∣
CÁLCULO / CIENCIAS AMBIENTALES / TEMA 5 93
hemos demostrado por tanto que el error puede acotarse por el máximo de la derivada cuarta de
la función.
Una consecuencia inmediata de este resultado es que si tenemos que integrar un
polinomio de grado 3, la integraciónexacta por la regla de Barrow y la “aproximada”
por el Método de Simpson (independientemente de la anchura del intervalo) coinciden,
el error es exactamente cero.
Una explicación gráfica de este sorprendente resultado (no olvidemos que Simpson se
basa en integrar un polinomio de grado 2, diferente por tanto al integrando, polinomio
de grado 3), la observamos en la Figura 6.1.
a xm b
Figura 5.1: Gráfica de un polinomio de grado 3 en un intervalo [a, b] y del correspondiente
polinomio de grado dos (en gris) que interpola los puntos de abscisa a, xm y b. Puede observarse
como el error de interpolación (por defecto) entre a y xm es idéntico al error (por exceso) entre
xm y b.
Método de Simpson Compuesto
De manera completamente análoga a lo expuesto para el Método de los Trapecios, es
posible generalizar (mejorando la precisión) el Método de Simpson por medio de la sub-
división del intervalo dado en otros más reducidos. De esta forma si partimos el intervalo
[a, b] en n subintervalos de anchura h = b−an tendremos la partición: {x0, x1, . . . , xn}. De
cara a aplicar el Método de Simpson simple paso a paso observamos inmediatamente que
n debe ser un número par para conseguir que todo [a, b] quede incluido en la integración
numérica. Tendremos entonces:∫ b
a
f(x) dx =
∫ x2
a
f(x)dx+
∫ x4
x2
f(x)dx+ . . .+
∫ xn
xn−2
f(x)dx
y los puntos x1, x3, . . . , xn−1 representarán el papel de “puntos medios” en cada una de
las aplicaciones sucesivas del método simple.
94 CÁLCULO / CIENCIAS AMBIENTALES / TEMA 5
De forma expĺıcita se obtiene:∫ b
a
f(x) dx ≈ h
3
(f(a) + 4I + 2P + f(b))
donde I y P representan las sumas:
I =
n−1∑
i=1, impares
f(xi) = f(x1) + f(x3) + . . .+ f(xn−1)
P =
n−2∑
i=2,pares
f(xi) = f(x2) + f(x4) + . . .+ f(xn−2)
De cara a la estimación del error, en cada uno de los pasos deberemos considerar
E ≤
∣∣∣∣h590M4
∣∣∣∣
De esta forma, el error de integración en el Método compuesto vendrá dado por:
E ≤
∣∣∣∣h590 (M14 +M24 + ...+M n24 )
∣∣∣∣ ≤ ∣∣∣∣h590 n2M4
∣∣∣∣
donde se denota M i4 a los máximos de la derivada cuarta en cada aplicación del método
simple y M4 al máximo de la derivada cuarta en todo [a, b].
Concluimos por tanto en la expresión:
E ≤
∣∣∣∣b− a180 h4M4
∣∣∣∣
Ejemplo 1. Calcular el valor aproximado de la integral∫ 1
0
x dx
(x+ 1)(x+ 2)
utilizando la regla de Simpson compuesta con n = 8.
Recordemos la tabla de valores utilizadas en la sección anterior al realizar este ejercicio
mediante el método de los trapecios:
x0 x1 x2 x3 x4 x5 x6 x7 x8
0. 0.125 0.25 0.375 0.5 0.625 0.75 0.875 1.0
f(x0) f(x1) f(x2) f(x3) f(x4) f(x5) f(x6) f(x7) f(x8)
0. 0.05228 0.08888 0.11483 0.13333 0.14652 0.15584 0.162319 0.16666
de manera que
I ≈ h
3
[f(x0) + 4 (f(x1) + f(x3) + f(x5) + f(x7)) + 2 (f(x2) + f(x4) + f(x6)) + f(x8)]
≈ 0.125
3
[4(0.05228 + 0.11483 + 0.14652 + 0.162319) + 2(0.0888 + 0.1333 + 0.15584) + 0.1666]
≈ 0.117773
CÁLCULO / CIENCIAS AMBIENTALES / TEMA 5 95
que al ser comparado con el valor exacto 0.1177830 y el obtenido por la regla de los trapecios
0.117166 nos permite concluir que este método es más preciso que el anterior.
Comparando de manera general los dos métodos tendremos:
n I(Trapecios) I(Simpson)
n = 1 0.08333
n = 2 0.108333 0.116667
n = 3 0.113492
n = 4 0.11535 0.117689
n = 5 0.11622
n = 6 0.117763
n = 8 0.117776
n = 10 0.11739 0.11778
n = 100 0.1177791 0.117783
Ejemplo 2. Teniendo en cuenta que no es conocida una primitiva de la función f(x) = ex
2
,
calcúlese el valor de la integral definida ∫ 1
0
ex
2
dx
con un error menor a 0.003.
La función con la que debemos trabajar es ex
2
. Aplicando la fórmula de Simpson cometemos
un error dado por
E ≤
∣∣∣∣(b− a) h4180M4
∣∣∣∣ , M4 ≥ f (4)(x) , ∀x ∈ [0, 1]
Calcularemos las derivadas correspondientes:
f ′(x) = 2xex
2
f ′′(x) = 2(1 + 2x2)ex
2
f ′′′(x) = 4(3x+ 2x3)ex
2
f (iv)(x) = 4(4x4 + 12x2 + 3)ex
2
Se puede observar que f (iv)(x) es creciente en [0, 1] de modo que el máximo valor de dicha función
coincide con el valor en x = 1, esto es, f (iv)(1) = 4e1(4 + 12 + 3) < 4 · 3 · 19 = 228, por lo que
consideraremos que M4 ≤ 228. Por ello,
E(N) ≤ (b− a)
5
180N4
228 =
19
15N4

E(1) ≈ 1.2666
E(2) ≈ 0.0791
E(3) ≈ 0.0156
E(4) ≈ 0.0049
E(5) ≈ 0.0020
de modo que para que el número de subintervalos sea par hemos de tomar
N = 6
96 CÁLCULO / CIENCIAS AMBIENTALES / TEMA 5
xi x0 x1 x2 x3 x4 x5 x6
= 0. 1/6 1/3 1/2 2/3 5/6 1.
f(xi) f(x0) f(x1) f(x2) f(x3) f(x4) f(x5) f(x6)
= 1. 1.02817 1.11752 1.28403 1.55962 2.0026 2.71828
Finalmente:
I ≈ h
3
[f(x0) + 4 (f(x1) + f(x3) + f(x5)) + 2 (f(x2) + f(x4)) + f(x8)]
≈ 1/6
3
[1 + 4(1.02817 + 1.28403 + 2.0026) + 2(1.11752 + 1.55962) + 2.71828]
≈ 1.4628
Ejemplo 3. Una cuerda vibra adoptando la forma,
y = senx
entre las abscisas x = 0 y x = 4 en un instante t0. Calcúlese aproximadamente la longitud de la
cuerda, utilizando un método numérico con n = 8.
Dado que tenemos que calcular la longitud de la función f(x) = senx, entre x = 0 y x = 4,
aplicaremos la fórmula
L =
∫ b
a
√
1 + (f ′(x))2dx =
∫ 4
0
√
1 + cos2 xdx
que nos proporciona la integral que debemos estimar numéricamente mediante la regla de Simpson
con n = 8 como propone el enunciado.
g(x) =
√
1 + cos2 x
xi x0 x1 x2 x3 x4 x5 x6 x7 x8
= 0. 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0
g(xi) g(x0) g(x1) g(x2) g(x3) g(x4) g(x5) g(x6) g(x7) g(x8)
= 1.41421 1.33047 1.13663 1.0025 1.08313 1.28134 1.40715 1.37002 1.19468
I ≈ h
3
[g(x0) + 4 (g(x1) + g(x3) + g(x5) + g(x7)) + 2 (g(x2) + g(x4) + g(x6)) + g(x8)]
≈ 0.5
3
[1.41421 + 4(1.33047 + 1.0025 + 1.28134 + 1.37002) + 2(1.13663 + 1.08313 + 1.40715) + 1.19468]
≈ 4.96667
Es posible calcular de forma precisa, por otros métodos, este resultado, obteniéndose: 4.966615,
por lo que deducimos que el Método de Simpson proporciona un valor muy correcto en este caso.
Ejemplo 4. Un agricultor desea conocer la superficie aproximada de un prado limitado por una
carretera, dos caminos perpendiculares a ella y la ribera de un ŕıo, de manera que si colocamos
CÁLCULO / CIENCIAS AMBIENTALES / TEMA 5 97
unos ejes cartesianos sobre la carretera (eje OX) y uno de los caminos (eje OY, abscisa x = 0),
el segundo camino será la recta vertical x = 2 (unidades en cientos de metros). Se toman varias
medidas desde la carretera hasta la ribera, obteniéndose las siguientes coordenadas para los
puntos de la ribera: (0, 1.5), (0.5, 1.8), (1, 2.1), (1.5, 1.75), (2, 1.3).
Calcular aproximadamente el área de dicho terreno utilizando las reglas de los trapecios y
de Simpson. Determinar el área si extendemos el terreno hasta la abscisa x = 2.5 sabiendo que
el ŕıo en tal caso pasa por el punto (2.5, 1.1).
En este caso desconocemos la función de forma expĺıcita, teniendo en cuenta tan solo los
valores de la tabla que nos han sido facilitados. Se tiene:
xi x0 x1 x2 x3 x4
= 0. 0.5 1.0 1.5 2.0
f(xi) f(x0) f(x1) f(x2) f(x3) f(x4)
= 1.5 1.8 2.1 1.75 1.3
de modo que usando el método de los trapecios podemos escribir
I1 ≈
h
2
[f(x0) + 2 (f(x1) + f(x2) + f(x3)) + f(x4)]
≈ 0.5
2
[1.5 + 2(1.8 + 2.1 + 1.75) + 1.3]
≈ 3.4333
mientras que si usamos el método de Simpson se llega a
I2 ≈
h
3
[f(x0) + 4 (f(x1) + f(x3)) + 2 (f(x2)) + f(x4)]
≈ 0.5
3
[1.5 + 4(1.8 + 1.75) + 2(2.1) + 1.3]
≈ 3.5333
Si se añade un nuevo punto, la tabla queda dada por
xi x0 x1 x2 x3 x4 x5
= 0. 0.5 1.0 1.5 2.0 2.5
f(xi) f(x0) f(x1) f(x2) f(x3) f(x4) f(x5)
= 1.5 1.8 2.1 1.75 1.3 1.1
Ahora la regla de los trapecios proporcionará:
I ≈ h
2
[f(x0) + 2 (f(x1) + f(x2) + f(x3) + f(x4)) + f(x5)]
≈ 0.5
2
[1.5 + 2(1.8 + 2.1 + 1.75 + 1.3) + 1.1] ≈ 4.125
mientras que si el método de Simpson no es aplicable de forma directa dado que estamos con-
siderando un número impar de subintervalos en este caso. Lo que podemos hacer es considerar
la regla de Simpson para los 4 subintervalos primeros y estimar el quinto subintervalo mediante
la regla de los trapecios. Aśı queda
I = I2 + I
′ = 3.5333+
0.5
2
(1.3 + 1.1) ≈ 4.13333