UNIVERSIDAD DE EXTREMADURA - matematicas.unex.esmatematicas.unex.es/~ghierro/Matematicas...
49
Transcript of UNIVERSIDAD DE EXTREMADURA - matematicas.unex.esmatematicas.unex.es/~ghierro/Matematicas...
UNIVERSIDAD DE EXTREMADURADepartamento de Matematicas
Matematicas
Manuel Fernandez Garcıa-Hierro
Badajoz, Febrero 2008
Capıtulo XIntegracion numerica
Introduccion
La integral definida
I(f) =∫ b
a
f(x) dx
se define como el lımite de las sumas de Riemann. Bajo ciertas
hipotesis sobre f(x) (por ejemplo, f(x) continua), se prueba
que
I(f) = F (b)− F (a)
donde F es una primitiva de f , es decir, es una funcion tal
que F ′(x) = f(x). Muchas integrales pueden ser calculadas
usando esta formula. Sin embargo, la mayor parte de las
integrales no pueden ser evaluadas mediante el procedimiento
anterior, porque la primitiva, F de f , no se puede expresar en
terminos de funciones elementales (polinomicas, exponenciales,
trigonometricas, sus inversas, composiciones de las funciones
anteriores,. . . ). Ejemplos de este tipo de integrales son∫ 1
0
e−x2dx,
∫ π
0
xπ sen(√x) dx.
Ası que se necesitan otros metodos para calcular estas
integrales.
En la primera seccion del capıtulo, se definen dos de los
metodos mas populares de calculo aproximado de I(f): la
regla de los trapecios y la de Simpson. En la segunda seccion
se analiza el error que se comete al usar dichos metodos.
La regla del trapecio y de Simpson
La regla del trapecio La idea central para calcular
aproximadamente
I(f) =∫ b
a
f(x) dx
es reemplazar f(x) por una funcion que aproxima a f(x) y
cuya integral sabemos calcular. En primer lugar vamos a usar
el polinomio lineal de interpolacion. De este modo aproximamos
f(x) por
P1(x) =(b− x)f(a) + (x− a)f(b)
b− a
que interpola a f(x) en a y b. Entonces
T1(f) = I(P1) =∫ b
a
P1(x) dx
= (b− a)(f(a) + f(b)
2
).
Si f(x) es parecida a un funcion lineal en [a, b], T1(f) deberıa
ser una buena aproximacion a I(f). T1(f) es el area del
trapecio sombreado que se muestra en la figura siguiente.
Ejemplo 1. Aproximar la integral
I =∫ 1
0
dx
1 + x.
El valor verdadero es I = ln 2 .= 0.693147. Usando la formula
del trapecio, se obtiene
T1 =12
(1 +
12
)=
34
= 0.75.
El error es
I − T1.= −0.0569
Para mejorar la aproximacion T1(f) cuando f(x) no es parecida
a una funcion lineal en [a, b], dividimos el intervalo [a, b] en
subintervalos mas pequenos a los que le aplicamos la regla del
trapecio. Si los subintervalos son suficientemente pequenos,
f(x) debe estar proxima a una funcion lineal en cada uno de
ellos. Ver la figura segunda.
Ejemplo 2. Calcular aproximadamente mediante la formula
del trapecio
I =∫ 1
0
dx
1 + x
usando T1(f) sobre dos subintervalos de igual longitud.
I =∫ 1/2
0
dx
1 + x+∫ 1
1/2
dx
1 + x
I.=
12
(1 + 2
3
2
)+
12
(23 + 1
2
2
)
T2 =1724
= 0.70833
I − T2.= −0.0152
Se observa que I − T2 es aproximadamente la cuarta parte de
I − T1
A continuacion se deduce una formula general que usa n
subintervalos de igual longitud. Sea
h =b− an
la longitud de cada subintervalo. Los extremos de los
subintervalos son
xj = a+ jh, j = 0, . . . , n
Entonces
I(f) =∫ b
a
f(x) dx =∫ xn
x0
f(x) dx
=∫ x1
x0
f(x) dx+∫ x2
x1
f(x) dx
+∫ xn
xn−1
f(x) dx.
Se aproxima cada integral mediante la regla del trapecio y se
tiene en cuenta que cada subintervalo [xj−1, xj] tiene longitud
h, para obtener
I(f) .= hf(x0) + f(x1)
2+ h
f(x1) + f(x2)2
+ hf(x2) + f(x3)
2+ · · ·+ h
f(xn−1) + f(xn)2
.
Y finalmente
Tn(f) =
h
(12f(x0) + f(x1) + · · ·+ f(xn−1) +
12f(xn)
)Esta formula se llama regla de integracion numerica del
trapecio. El subındice n indica el numero de subintervalos
que se utilizan. Los puntos x0, . . . , xn se llaman nodos de
integracion numerica.
Antes de presentar algunos ejemplos de Tn(f) discutiremos
sobre la eleccion de n. Usualmente una sucesion creciente de
valores de n hara que Tn(f) se aproxime cada vez mas a I(f).
Si elegimos la sucesion
T2(f), T4(f), T8(f), . . .
resulta que el calculo de T2n(f) utiliza todos los valores de
f(x) del computo anterior. Para ilustrar la manera en la que
los valores son reutilizados cuando se duplica n, consideremos
el calculo de T2(f) y T4(f).
T2(f) = h
(f(x0)
2+ f(x1) +
f(x2)2
)con
h =b− a
2, x0 = a
x1 =a+ b
2x2 = b.
Tambien
T4(f) = h
(f(x0)
2+ f(x1) + f(x2) + f(x3) +
f(x4)2
)
con
h =b− a
4x0 = a
x1 =3a+ b
4x2 =
a+ b
2
x3 =a+ 3b
4x4 = b.
En la formula para T4(f), solamente f(x1) y f(x3) necesitan
ser calculados, porque el resto de valores ya es conocido a partir
del calculo de T2(f). Por esta razon en todos los ejemplos que
siguen se duplica n.
Ejemplo 3. Calculamos las integrales
I(1) =∫ 1
0
e−x2dx
.= 0.74682413281234
I(2) =∫ 4
0
dx
1 + x2= tan−1(4) .= 1.3258176636680
I(3) =∫ 2π
0
dx
2 + cosx=
2π√3.= 3.6275987284684
Los resultados se muestran en la tabla que sigue. Solo se dan
los errores I(f)−Tn(f), puesto que es la magnitud que interesa
conocer si queremos determinar la velocidad de convergencia de
Tn(f) a I(f). La columna con la etiqueta ”Ratio” da la razon
entre errores sucesivos, es decir, el factor de decrecimiento del
error cuando n se duplica.
De la tabla se observa que el error al calcular I(1) e I(2) decrece
con un factor de 4 cuando n se duplica. En el tercer ejemplo
se observa que I(3) decrece mucho mas rapidamente.
La regla de Simpson Con el fin de mejorar el error
cometido al calcular T1(f) en lugar de I(f), se utiliza el
polinomio cuadratico de interpolacion para aproximar f(x) en
el intervalo [a, b].
Sea P2(f) el polinomio cuadratico que interpola a f(x) en
a, c = (a+ b)/2 y b. Se obtiene
I(f) .=∫ b
a
P2(x) dx
=∫ b
a
((x− c)(x− b)(a− c)(a− b)
f(a) +(x− a)(x− b)(c− a)(c− b)
f(c)
+(x− a)(x− c)(b− a)(b− c)
f(b))dx
Esta integral se puede calcular directamente, pero es mas facil
llamar
h =b− a
2
y hacer el cambio de variable u = x−a en la integral. Entonces
∫ b
a
(x− c)(x− b)(a− c)(a− b)
dx
=1
2h2
∫ a+2h
a
(x− c)(x− b) dx
=1
2h2
∫ 2h
0
(u− h)(u− 2h) du
=1
2h2
[u3
3− 3
2u2h+ 2h2u
]2h0
=h
3.
De donde se obtiene
S2(f) =h
3
(f(a) + 4f
(a+ b
2
)+ f(b)
)
El metodo se ilustra en la figura siguiente
Ejemplo 4. Utilice la integral
I =∫ 1
0
dx
1 + x
para calcular la aproximacion S2. En este caso h = (b−a)/2 =12, y
S2 =1/23
(1 + 4
(23
)+
12
)=
2536
.= 0.69444.
El error es
I − S2 = ln 2− S2.= −0.00130.
Para comparar este error con el error en la regla de los trapecios,
se debe usar T2, puesto que el numero de evaluaciones sobre
f(x) es el mismo tanto para S2 como para T2. Puesto que
I − T2.= −0.0152
el error usando Simpson es aproximadamente 12 veces menor
que el error cometido usando la regla de los trapecios.
La regla de Simpson S2(f) sera una aproximacion precisa
de I(f) si f(x) es parecida a una funcion cuadratica en [a, b].
Cuando no es ası, procedemos de la misma forma que en la regla
de los trapecios. Sea n un entero positivo par, h = (b − a)/ny definimos los puntos de evaluacion para f(x) por
xj = a+ jh, j = 0, 1, . . . , n.
Dividimos [a, b] en los subintervalos
[x0, x2], [x2, x4], . . . [x2n−2, x2n].
Notese que cada uno de ellos contiene tres puntos de evaluacion
para f(x). Entonces
I(f) =∫ b
a
f(x) dx =∫ x2n
x0
f(x) dx
=∫ x2
x0
f(x) dx+∫ x4
x2
f(x) dx
+ · · ·+∫ x2n
x2n−2
f(x) dx.
Aproximando cada integral por Simpson
I(f) .=h
3(f(x0) + 4f(x1) + f(x2))
+h
3(f(x2) + 4f(x3) + f(x4))
+ · · ·+ h
3(f(xn−2) + 4f(xn−1) + f(xn)) .
Simplificando
Sn(f) =h
3(f(x0) + 4f(x1)
+ 2f(x2) + 4f(x3) + 2f(x4)
+ · · ·+ 2f(xn−2 + 4f(xn−1) + f(xn))).
A la formula anterior la llamaremos regla de Simpson. El ındice
n es el numero de subdivisiones utilizadas para definir los nodos
de integracion x0, x1, . . . , xn.
Ejemplo 5. Calculamos las integrales
I(1) =∫ 1
0
e−x2dx
.= 0.74682413281234
I(2) =∫ 4
0
dx
1 + x2= tan−1(4) .= 1.3258176636680
I(3) =∫ 2π
0
dx
2 + cosx=
2π√3.= 3.6275987284684
Los resultados usando la regla de Simpson se muestran en la
tabla siguiente
Se observa que para las integrales I(1) e I(2) el factor de
disminucion del error es de aproximadamente 16, mientras que
I(3) converge mucho mas rapidamente.
Formulas del error
Formula del error para la regla del trapecio
En todo este apartado supondremos que las funciones son dos
veces derivables con continuidad en su intervalo de definicion.
Se cumple que
ETn (f) := I(f)− Tn(f) = −h2(b− a)
12f ′′(cn),
donde h = (b− a)/n y cn ∈ [a, b].
la formula anterior se utiliza para acotar el error ETn (f),
sustituyendo f′′(cn) por su mayor valor posible en el intervalo
[a, b]. Ası se hara en el ejemplo que sigue. Notese tambien que
la formula del error es consistente con los errores observados
para las integrales I(1) e I(2) en la Tabla 1. Cuando n se
duplica, h se divide por 2 y h2 se divide por 4. Este es el
factor de decrecimiento para el error trapezoidal observado en
la Tabla 1.
Ejemplo 6. Sea
I =∫ 1
0
dx
1 + x= ln 2.
Aquı
f(x) =1
1 + x, f ′′(x) =
2(1 + x3)
.
Sustituyendo en la formula del error trapezoidal, se obtiene
ETn (f) = −h2
12f ′′(cn), 0 ≤ cn ≤ 1, h =
1n.
Esta formula no puede ser calculada de forma exacta, porque
cn es desconocido. Pero podemos acotar el error de la siguiente
manera: Puesto que
max0≤x≤1
2(1 + x)3
= 2,
se tiene
|ETn (f)| ≤ h2
12(2) =
h2
6.
Para n = 1 y n = 2, se obtiene
|ET1 (f)| ≤ 16.= 0.167, |ET2 (f)| ≤ (1/2)2
6.= 0.0417.
Comparando estas cotas con los errores verdaderos
I − T1.= −0.0569, I − T2
.= −0.0152,
se observa que son de dos a tres veces mayores.
Una estimacion del error de la formula del
trapecio La demostracion de la formula del error
ETn (f) := I(f)− Tn(f) = −h2(b− a)
12f ′′(cn),
se fundamenta en
∫ α+h
α
f(x) dx− h(f(α) + f(α+ h)
2
)(1)
= −h3
12f ′′(c), (2)
donde c ∈ [α, α + h]. Supuesta probada esta formula,
demostraremos la formula del error.
ETn (f) =∫ b
a
f(x) dx− Tn(f)
=∫ xn
x0
f(x) dx− Tn(f)
=∫ x1
x0
f(x) dx− h(f(x0) + f(x1)
2
)+∫ x2
x1
f(x) dx− h(f(x1) + f(x2)
2
)+ . . .
+∫ xn
xn−1
f(x) dx− h(f(xn−1) + f(xn)
2
)
Aplicando 1 a cada uno de los sumandos anteriores, se obtiene
ETn (f)
= −h3
12f ′′(t1)−
h3
12f ′′(t2)− · · · −
h3
12f ′′(tn)
= −h2
12(hf ′′(t1) + hf ′′(t2) + · · ·+ hf ′′(tn))
= −h2
12(b− a)f
′′(t1) + f ′′(t2) + · · ·+ f ′′(tn)n
= −h2
12(b− a)f ′′(cn)
Las constantes desconocidas t1, . . . , tn estan localizadas en los
intervalos
[x0, x1], [x1, x2], . . . , [xn−1, xn]
respectivamente y cn ∈ [a, b].
La ultima igualdad es consecuencia del Teorema de los
valores intermedios para funciones continuas, cuyo enunciado
presentamos a continuacion.
Teorema. Sea f(x) una funcion continua en el intervalo [a, b].Sean
m = minx∈[a,b]
f(x), M = maxx∈[a,b]
f(x).
Entonces para cada valor y0 que cumple m ≤ y0 ≤ M , existe
un valor c ∈ [a, b] tal que f(c) = y0.
Ahora todo lo que hay que hacer es aplicar este resultado a la
funcion f ′′(x) que suponemos continua en [a, b]. En efecto,
m = minx∈[a,b]
f ′′(x),
M = maxx∈[a,b]
f ′′(x).
Entonces
m ≤ f ′′(t1) + f ′′(t2) + · · ·+ f ′′(tn)n
≤M.
Por tanto existe cn ∈ [a, b] tal que
f ′′(cn) =f ′′(t1) + f ′′(t2) + · · ·+ f ′′(tn)
n.
Para obtener una estimacion del error trapezoidal, observese
que
hf ′′(t1) + hf ′′(t2) + · · ·+ hf ′′(tn)es una suma de Riemann para la integral∫ b
a
f ′′(x) dx = f ′(b)− f ′(a)
y que cuando n→∞ estas sumas de Riemann convergen a la
integral. Entonces
ETn (f)
= −h2
12(hf ′′(t1) + hf ′′(t2) + · · ·+ hf ′′(tn))
.= −h2
12
∫ b
a
f ′′(x) dx
= −h2
12(f ′(b)− f ′(a)) = ETn (f).
Esta estimacion del error se dice asintotica ya que mejora
conforme n aumenta. Es tan facil de calcular como lo sea f ′.
Ejemplo 7. Sea
I =∫ 1
0
dx
1 + x= ln 2.
Entonces f ′(x) = −1/(1 + x)2, h = 1/n y
ETn (f) = −h2
12
(−1
(1 + 1)2− −1
(1 + 0)2
)=−h2
16.
Para n = 1 y n = 2 se tiene
ET1 (f) = − 116
= −0.0625, ET2 (f) = −0.0156,
que son bastante proximos a los errores
I − T1.= −0.0569, I − T2
.= −0.0152.
La estimacion ETn (f) tiene varias ventajas sobre ETn (f). En
primer lugar confirma que cuando n se duplica, el error
disminuye en un factor de 4, supuesto que f ′(b) − f ′(a) 6= 0.
Esto coincide con los resultados para I(1) e I(2) de la Tabla 1.
En segundo lugar, se observa que la convergencia de Tn(f) es
mas rapida cuando f ′(b) − f ′(a) = 0. Esto explica la mayor
velocidad de convergencia para I(3) en la Tabla 1. Finalmente,
la estimacion ETn (f) conduce a una formula de integracion mas
precisa. Para ello basta tener en cuenta:
I(f)− Tn(f) .=−h2
12(f ′(b)− f ′(a))
I(f) .= Tn(f)− h2
12(f ′(b)− f ′(a)) = CTn(f).
Esta ultima formula se llama regla del trapecio corregida.
Formula del error para la regla de Simpson El
tipo de analisis utilizado en la seccion anterior tambien sirve
para calcular el error y el error asintotico al usar la regla de
Simpson.
Supongase que f(x) es cuatro veces derivable con
continuidad en [a, b]. Entonces
ESn(f) = I(f)− Sn(f) = −h4(b− a)180
f (4)(cn),
donde cn ∈ [a, b] y h = (b − a)/n. Ademas, este error puede
estimarse con la formula asintotica del error
ESn(f) .= − h4
180(f ′′′(b)− f ′′′(a)).
Notese que la primera formula dice que la regla de Simpson es
exacta para todas las f(x) que son polinomios de grado ≤ 3,
mientras que la interpolacion cuadratica en la que esta basada
la regla de Simpson es exacta solo cuando f(x) es un polinomio
de grado ≤ 2. este hecho, junto a la alta potencia de h en
la formula del error y la sencillez del metodo, hacen que la
formula de Simpson sea el metodo mas popular de integracion
numerica.
Al igual que se hizo con la regla del trapecio, se define la
regla de Simpson corregida de la manera siguiente.
CSn(f) = Sn(f)− h4
180(f ′′′(b)− f ′′′(a)).
REFERENCIAS
Kendall Atkinson, Elementary Numerical Analysis, Second
Edition, Jhon Wiley & Sons, Inc. New York, 1993
A. Cordero Barbero, J.L. Hueso Pagoaga, E. Martınez
Molada, J.R. Torregrosa Sanchez, Problemas resueltos de
Metodos Numericos, Edit. Thomson, Madrid, 2006, ISBN
84-9732-409-9.
