Cap tulo IV Integraci on Num erica - Sede Medellinearamosn/cursos/... · 2013. 10. 28. · Cap tulo...

Capıtulo IV

Integracion Numerica

IV.1. Cuadraturas: Reglas Simples

La formulas de cuadratura o reglas simples se obtienen por medio de interpo-lacion polinomial: la funcion a integrar se “muestrea”, es decir, se toman puntosde la funcion

(x0, y0), (x1, y1), . . . , (xn, yn)

donde yk = f(xk) para k = 0, . . . , n, se determina el polinomio interpolante pn(x)de grado a lo mas n para esos puntos, y finalmente se integra. Intuitivamente, sipn(x) es una buena aproximacion de f(x), entonces su integral

∫ b

apn(x)dx debe

ser una buena aproximacion de∫ b

af(x)dx.

Regla del Trapecio

Para comenzar, solo dos puntos (a, f(a)) y (b, f(b)) e interpolacion lineal re-sulta en ∫ b

a

f(x)dx ≈ b− a2· (f(a) + f(b))

Esta es la regla del trapecio: aproxima la integral por medio del area bajo eltrapecio determinado por la interpolacion lineal.

Caso General

Para el caso general, es conveniente usar el metodo de interpolacion de La-grange. Para los puntos (x0, y0), . . . , (xn, yn), con yk = f(xk) el polinomio deinterpolacion esta dado por

pn(x) =n∑

k=0

yk · Ln,k(x)

donde

Ln,k(x) =

∏j 6=k(x− xj)∏j 6=k(xk − xj)

IV.1

IV.2 CAPITULO IV. INTEGRACION NUMERICA

De esta representacion, integrando, obtenemos∫ b

a

f(x)dx ≈∫ b

a

pn(x)dx

=

∫ b

a

(n∑

k=1

yk · Ln,k(x)

)dx

=n∑

k=0

(∫ b

a

Ln,k(x)dx

)yk

=n∑

k=0

wkyk

donde

wk =

∫ b

a

Ln,k(x)dx.

Note que estos pesos wk dependen solo de los nodos, no de la funcion. La suma

Q[f ] =n∑

k=0

wkfk

se llama una formula de cuadratura.

Regla de Simpson

Consideremos el caso de interpolacion cuadratica. Se tienen tres nodos x0 = a,x1 = (a + b)/2 y x2 = b. Por conveniencia, sea h = (b − a)/2, de tal forma quexk = a+ kh, k = 0, 1, 2.

Entonces los tres pesos w0, w1, w2 estan dados por

w0 =∫ x2

x0

(x− x1)(x− x2)(x0 − x1)(x0 − x2)

dx =∫ 2h

0

(x− h)(x− 2h)(−h)(−2h)

dx

=1

2h2

[x3

3− x2

23h + 2h2x

]2h

0

=13h

w1 =∫ x2

x0

(x− x0)(x− x2)(x1 − x0)(x1 − x2)

dx =∫ 2h

0

x(x− 2h)(h)(−h)

dx = − 1h2

[x3

3− x2h

]2h

0

=43h

w2 =∫ x2

x0

(x− x0)(x− x1)(x2 − x0)(x2 − x1)

dx =∫ 2h

0

x(x− h)(2h)(h)

dx =1

2h2

[x3

3− x2

2h

]2h

0

=13h

(las segundas integrales se han obtenido por medio de una translacion). Unaforma alternativa de obtenerlos wk, que resulta mas facil de sistematizar paraobtener reglas de mas alto orden, es la siguiente. Primero, como la translacionde la integral arriba lo demuestra, podemos asumir que el intervalo de intereses [a, b] = [0, 2h] y xk = kh, k = 0, 1, 2. Como estamos usando interpolacion

IV.1. CUADRATURAS: REGLAS SIMPLES IV.3

cuadratica, la interpolacion y por lo tanto la integral obtenida con la formula decuadratura deben ser exactas. Ası que para f(x) = xk, k = 0, 1, 2 se debe tenerque, para j = 0, 1, 2,

2∑k=0

wkxjk =

∫ 2h

0

xjdx =

[xj+1

j + 1

]2h

0

=(2h)j+1

j + 1

Estas tres ecuaciones explıcitamente son

w0 +w1 +w2 = 2h+w1h +w2(2h) = 2h2

+w1h2 +w2(2h)2 = 8

3h3

y cancelando h’sw0 +w1 +w2 = 2h

+w1 +2w2 = 2h+w1 +4w2 = 8

3h

La solucion es w0 = h/3, w1 = 4h/3 y w2 = h/3, igual a lo obtenido por mediode las integrales anteriormente. La cuadratura que se obtiene se llama la regla deSimpson: ∫ x2

x0

f(x)dx ≈ h

3(f0 + 4f1 + f2)

Solucion General

Generalizando, con xk = kh, k = 0, 1, . . . , n, obtenemos, para j = 0, 1, . . . , n,

n∑k=0

wkxjk =

∫ nh

0

xjdx =

[xj+1

j + 1

]nh

0

=(nh)j+1

j + 1

y de aquı,n∑

k=0

wk(kh)j =(nh)j+1

j + 1.

Cancelando hj y con wk = wk/h, se obtiene el sitema de ecuaciones

n∑k=0

wkkj =

nj+1

j + 1, para j = 0, 1, . . . , n,

que en forma expandida es1 1 1 · · · 10 1 2 · · · n0 1 4 · · · n2

......

.... . .

...0 1 2n · · · nn

w0

w1

w2...wn

=

nn2

2n3

3...

nn+1

n+1


Note que la matriz de coeficientes es una matriz de Vandermonde. Podemos usarMatlab para solucionar estos sistemas de ecuaciones. Con la siguiente funcion

function cuadratura(N)for n=1:N

nn=(n+1):-1:1;b=((n.^nn)./nn)’;V= (vander(0:1:n))’;c= (V\b)’

end

obtenemos

>> cuadratura(6)c =

1/2 1/2c =

1/3 4/3 1/3c =

3/8 9/8 9/8 3/8c =

14/45 64/45 8/15 64/45 14/45c =

95/288 125/96 125/144 125/144 125/96 95/288c =

41/140 54/35 27/140 68/35 27/140 54/35 41/140

Las dos primeras soluciones son las reglas del trapecio y de Simpson. Las dossiguientes se llaman reglas de Simpson de 3/8 y de Boole:

Regla de Simpson-3/8

∫ x3

x0

f(x)dx =3h

8(f0 + 3f1 + 3f2 + f3)

Regla de Boole

∫ x4

x0

f(x)dx =2h

45(7f0 + 32f1 + 12f2 + 32f3 + 7f4)

Ejemplo.

En este ejemplo deseamos aproximar ln 5 por medio de la integral∫ 5

1dxx

. Lasgraficas muestran las interpolaciones de grados 1 a 4. Los resultados obtenidoscon Matlab muestran para cada valor de n, el valor de h y de la cuadratura, y unatabla con los valores del ındice k y xk, yk wk y wkyk. El valor de la cuadratura es,para cada n, la suma de la ultima columna (valores wkyk) multiplicada por h.

IV.1. CUADRATURAS: REGLAS SIMPLES IV.5

Figura IV.1: Integracion con las reglas del trapecio, Simpson, Simpson 3/8 y Boole.

>> cuadraturaf(8)n = 1h = 4cuad = 12/5[k x y w w*y] =

0 1 1 1/2 1/21 5 1/5 1/2 1/10

n = 2h = 2cuad = 76/45[k x y w w*y] =

0 1 1 1/3 1/31 3 1/3 4/3 4/92 5 1/5 1/3 1/15

n = 3h = 4/3cuad = 636/385[k x y w w*y] =

0 1 1 3/8 3/81 7/3 3/7 9/8 27/562 11/3 3/11 9/8 27/883 5 1/5 3/8 3/40

n = 4h = 1cuad = 364/225[k x y w w*y] =


0 1 1 14/45 14/451 2 1/2 64/45 32/452 3 1/3 8/15 8/453 4 1/4 64/45 16/454 5 1/5 14/45 14/225

n = 5h = 4/5cuad = 2083/1290[k x y w w*y] =

0 1 1 95/288 95/2881 9/5 5/9 125/96 625/8642 13/5 5/13 125/144 625/18723 17/5 5/17 125/144 625/24484 21/5 5/21 125/96 625/20165 5 1/5 95/288 19/288

n = 6h = 2/3cuad = 662/411[k x y w w*y] =

0 1 1 41/140 41/1401 5/3 3/5 54/35 162/1752 7/3 3/7 27/140 81/9803 3 1/3 68/35 68/1054 11/3 3/11 27/140 81/15405 13/3 3/13 54/35 162/4556 5 1/5 41/140 41/700

n = 7h = 4/7cuad = 2752/1709[k x y w w*y] =

0 1 1 1073/3527 1073/35271 11/7 7/11 810/559 509/5522 15/7 7/15 343/640 600/23993 19/7 7/19 649/536 1113/24954 23/7 7/23 649/536 646/17535 27/7 7/27 343/640 467/33616 31/7 7/31 810/559 551/16847 5 1/5 1073/3527 193/3172

n = 8h = 1/2cuad = 433/269[k x y w w*y] =

0 1 1 499/1788 499/17881 3/2 2/3 1183/712 1183/10682 2 1/2 -353/1348 -353/26963 5/2 2/5 388/131 776/6554 3 1/3 -1317/1028 -439/10285 7/2 2/7 388/131 765/9046 4 1/4 -353/1348 -353/53927 9/2 2/9 1183/712 1183/32048 5 1/5 499/1788 297/5321

Finalmente, la tabla siguiente muestra las cuadraturas y el error resultantes paran = 1, 2, . . . , 12.

IV.2. ERROR DE LAS REGLAS SIMPLES IV.7

n cuadratura error

1 2.400000000000000 0.7905620875659002 1.688888888888889 0.0794509764547883 1.651948051948052 0.0425101395139524 1.617777777777775 0.0083398653436755 1.614728626493345 0.0052907140592456 1.610704216419929 0.0012663039858297 1.610298283118095 0.0008603706839958 1.609665126452762 0.0002272140186619 1.609598154738391 0.00016024230429110 1.609482778174419 0.00004486574031811 1.609467583016121 0.00002967058202012 1.609424211461854 0.000013700972246

IV.2. Error de las Reglas Simples

Dada una cuadratura

Q[f ] =M∑

k=0

wkf(xk),

el error de truncamiento es

E[f ] =

∫ b

a

f(x)dx−Q[f ]

IV.2.1. Grado de Precision

El grado de precision de una formula de cuadratura es el entero positivo n talque E[p] = 0 para todo polinomio p(x) de grado ≤ n pero E[pn+1] 6= 0 para algunpolinomio pn+1(x) de grado n + 1. Veamos el grado de precision de las reglassimples antes derivadas. En cada caso, obviamente, el grado de precision es almenos el grado de la interpolacion para el cual se derivo la regla.

Trapecio

Por una parte ∫ h

0

x2dx =h3

3.

Por otra parte

QT [x2] =h

2(0 + h2) =

h3

2

Puesto que estos dos resultados son diferentes, el grado de precision de la regladel trapecio es 1.


Simpson

Tenemos que∫ 2h

0x3dx = (2h)4

4= 4h4 QS[x3] = h

3(0 + 4h3 + (2h)3) = 4h4∫ 2h

0x4dx = (2h)5

5= 32

5h5 QS[x4] = h

3(0 + 4h4 + (2h)4) = 20

3h5

Ası que el grado de la regla de Simpson es 3. Note que esto es uno mas de lo quepodıamos esperar dado que se obtuvo con interpolacion cuadratica.

Simpson 3/8

Tenemos que∫ 3h

0x4dx = (3h)5

5= 243

5h5; QS3/8[x

4] = 3h8

(0 + 3h4 + 3(2h)4 + (3h)4) = 992h5

Ası que el grado de la regla de Simpson 3/8 es 3.

Boole

Tenemos que∫ 4h

0x5dx = (4h)6

6= 2048

3h6; QB[x5] = 2h

45(7(0) + 32h5 + 12(2h)5 + 32(3h)5 + 7(4h)5)

= 20483h6∫ 4h

0x6dx = (4h)7

7= 16384

7h7; QB[x6] = 2h

45(7(0) + 32h6 + 12(2h)6 + 32(3h)6 + 7(4h)6)

= 70403h7

Ası que el grado de la regla de Boole es 5. De nuevo uno mas de lo que se podrıaesperar dado que el metodo se deriva de interpolacion de grado 4.

IV.2.2. Reglas Simples con Error

Se puede usar el error de interpolacion para estimar el error de las reglas decuadratura anteriores. Sin embargo, procediendo de esta manera no se obtiene lamejor estimacion posible. Por ejemplo, para la regla de Simpson se obtendrıa queel error es proporcional a f (3)(c)h4 para algun c ∈ [x0, x2]. Pero hemos obtenidoque esta regla es exacta para polinomios de grado 3, y consistente con esto se puedederivar (aunque no lo hacemos aquı) que el error es proporcional a f (4)(c)h5 paraalgun c ∈ [x0, x2], como se enuncia en esta seccion.

Regla del Trapecio

∫ x1

x0

f(x)dx =h

2(f0 + f1)−

h3

12f (2)(c) para algun c ∈ [x0, x1]

IV.2. ERROR DE LAS REGLAS SIMPLES IV.9

Prueba: El error de interpolacion lineal en x esta dado por

E1(x) = f(x)− p1(x) =1

2f ′′(c(x))(x− x0)(x− x1)

para algun c(x) ∈ [x0, x1]. Entonces

ET [f ] =

∫ x1

x0

1

2f ′′(c(x))(x− x0)(x− x1)dx

=1

2f ′′(c)

∫ x1

x0

(x− x0)(x− x1)dx

=1

2f ′′(c)

∫ h

0

x(x− h)dx

=1

2f ′′(c)

[1

3x3 − 1

2x2h

]h

0

= − 1

12f ′′(c)h3

donde en la segunda lınea se ha usado el teorema del valor medio para integralescon una funcion de peso (la funcion f ′′(c(x)) puede salir de la integral como f ′′(c)donde c = c(ξ) para algun ξ ∈ [x0, x1]).

Regla de Simpson∫ x2

x0

f(x)dx =h

3(f0 + 4f1 + f2)−

h5


Regla de Simpson-3/8∫ x3

x0

f(x)dx =3h

8(f0 + 3f1 + 3f2 + f3)−

3h5


Regla de Boole∫ x4

x0

f(x)dx =2h

45(7f0+32f1+12f2+32f3+7f4)−

8h7


Ejemplo.

Consideramos de nuevo la aproximacion de ln 5 por medio de la integral∫ 5

1dxx

.Se quiere determinar las cotas de error dadas teoricas de acuerdo a las formulasanteriores. Tenemos

f(x) =1

x; f ′(x) = − 1

x2 ; f ′′(x) = 2x3 ; f (3)(x) = − 6

x4

f (4)(x) = 24x5 ; f (5)(x) = −120

x6 ; f (6)(x) = 720x7


Todas estas son decrecientes en valor absoluto y por lo tanto el maximo ocurre enx = 1, en cada caso. Ası que

|ET (1/x, [1, 5])| ≤ h3

12maxx∈[1,5]

|f (2)(x)| = 43

12· 2 =

32

3= 10,666666

|ES(1/x, [1, 5])| ≤ h5

90maxx∈[1,5]

|f (4)(x)| = 25

90· 24 =

384

45= 8,533333

|ES3/8(1/x, [1, 5])| ≤ 3h5

80maxx∈[1,5]

|f (4)(x)| = 3(4/3)5

80· 24 =

512

135= 3,792592

|EB(1/x, [1, 5])| ≤ 8h7

945maxx∈[1,5]

|f (6)(x)| = 8

945· 720 =

128

21= 6,095238

Por otra parte, el valor exacto de la integral es ln 5 = 1,609437912434100 · · · yel resultado de las reglas del trapecio, de Simpson y de Boole son: 2,4, 1,6888888 · · ·y 1,61777777 · · ·. Ası que se observa que la cota de error esta lejos del error real.Esto se debe principalmente a que las derivadas decrecen rapidamente para x > 1pero tomamos el maximo en x = 1 para la cota.

IV.3. Reglas Simples Abiertas

Las reglas anteriores son cerradas porque incluyen los extremos del intervalocomo nodos. En las reglas abiertas se omiten los extremos. Las siguientes son lasprimeras tales reglas, para n nodos, con n = 1, 2, 3, 4, donde en cada caso a = x−1,b = xn, y x−1, x0, x1, . . . , xn dividen [a, b] en intervalos iguales:∫ x1

x−1

f(x)dx = 2hf(x0) +1

3f ′′(c)h3∫ x2

x−1

f(x)dx =3h

2(f(x0) + f(x1)) +

3

4f ′′(c)h3∫ x3

x−1

f(x)dx =4h

3(2f(x0)− f(x1) + 2f(x2)) +

14

45f (4)(c)∫ x4

x−1

f(x)dx =5h

24(11f(x0) + f(x1) + f(x2) + 11f(x0)) +

95

144f (4)(c)

IV.4. Reglas Compuestas

Como ocurre en interpolacion, en general se prefiere dividir el intervalo deintegracion en un numero de subintervalos y aplicar una regla simple de bajo ordenen cada uno de ellos. Las formulas resultantes se denominan reglas compuestas. Sea[a, b] el intervalo de integracion y M el numero de subintervalos. Para las reglasdel trapecio, de Simpson y de Boole se obtienen las siguientes reglas compuestas.Escribimos fk = f(xk).

IV.4. REGLAS COMPUESTAS IV.11

Trapecio

Sea h = (b − a)/M y xk = a + kh, k = 0, 1, . . . ,M . Se aplica la regla simpledel trapecio en cada subintervalo [xk, xk+1] y se obtiene la regla compuesta

T (f, h) =h

2(f0 + 2f1 + 2f2 + 2f3 + · · ·+ 2fM−2 + 2fM−1 + fM)

=h

2

(f0 +

M−1∑k=1

fk + fM

)

ET (f, h) = −(b− a)f (2)(c)

12h2 = O(h2) para algun c ∈ [a, b]

Simpson (1/3)

Sea h = (b− a)/2M y xk = a+ kh, k = 0, 1, . . . , 2M . Se aplica la regla simpledel trapecio en cada subintervalo [x2k, x2k+2] y se obtiene la regla compuesta

S(f, h) =h

3((f0 + 4f1 + f2) + (f2 + 4f3 + f4) + · · ·+ (f2M−2 + 4f2M−1 + f2M))

=h

3(f0 + 4f1 + 2f2 + 4f3 + 2f4 + · · ·+ 2f2M−2 + 4f2M−1 + f2M)

=h

3

(f(x0) + 2

M−1∑k=1

f2k + 4M−1∑k=0

f2k+1 + f2M

)

ES(f, h) = −(b− a)f (4)(c)

180h4 = O(h4) para algun c ∈ [a, b]

Boole

Sea h = (b− a)/4M y xk = a+ kh, k = 0, 1, . . . , 4M . Se aplica la regla simpledel trapecio en cada subintervalo [x4k, x4k+4] y se obtiene la regla compuesta

B(f, h) =2h

45

M−1∑k=0

(7f4k + 32f4k+1 + 12f4k+2 + 32f4k+3 + 7f4k+4)

=2h

45

(7f0 + 14

M−1∑k=1

f4k + 32M−1∑k=0

(f4k+1 + f4k+3) + 12M−1∑k=0

f4k+2 + 7f4M

)EB(f, h) = O(h6)

Ejemplo.

Se quiere aproximar la integral∫ 5

0exp(−x2/2)dx por medio de (i) la regla

compuesta del trapecio con M = 10 y de (ii) la regla compuesta de Simpson con


2M = 10. Para cada caso, (iii) cual es la cota teorica de error ? (iv) cual debe serel mınimo valor de M para garantizar error no mayor que 5× 10−9 ?

Sea f(x) = exp(−x2/2).

(i) Con M = 10, tenemos xk = k/2, para k = 0, 1, . . . , 10. Tenemos la tabla devalores

k xk fk

0 0 1,000000001 0,5 0,882496902 1,0 0,606530653 1,5 0,324652464 2,0 0,135335285 2,5 0,043936936 3,0 0,011108997 3,5 0,002187498 4,0 0,000335469 4,5 0,00004006

10 5,0 0,00000372

Entonces

T (f, 0,5) =0,5

2(1 + 2 · 0,88249690 + 2 · 0,60653065 + 2 · 0,32465246

+2 · 0,13533528 + 2 · 0,04393693 + 2 · 0,01110899 + 2 · 0,00218749

+2 · 0,00033546 + 2 · 0,00004006 + 0,00000372)

= 0,25 · 5,01325225 = 1,25331306

(ii) Usando la misma tabla de valores, tenemos

S(f, 0,5) =0,5

3(1 + 4 · 0,88249690 + 2 · 0,60653065 + 4 · 0,32465246

+2 · 0,13533528 + 4 · 0,04393693 + 2 · 0,01110899 + 4 · 0,00218749

+2 · 0,00033546 + 4 · 0,00004006 + 0,00000372)

= 0,16666667 · 7,51987997 = 1,25331333

(iii) Para usar las cotas teoricas, obtenemos las derivadas de f(x):

f (0)(x) = exp(−x2/2)

f (1)(x) = (−x) exp(−x2/2)

f (2)(x) = (−1 + x2) exp(−x2/2)

f (3)(x) = (3x− x3) exp(−x2/2)

f (4)(x) = (3− 6x2 + x4) exp(−x2/2)

f (5)(x) = (−15x+ 10x3 − x5) exp(−x2/2)

IV.4. REGLAS COMPUESTAS IV.13

Para la regla del trapecio

|ET (f, h)| ≤ (b− a)h2

12maxx∈[0,5]

|f (2)(x)|

Para determinar el max de |f (2)(x)| en [1, 5], encontramos los puntos crıticosde f (2)(x) igualando a cero f (3)(x):

f (3)(x) = (3x− x3) exp(−x2/2) = 0 ⇒ x = 0 o x =√

3.

Comparando entonces f(0), f(√

3) y f(5), se encuentra que el maximo (envalor absoluto) esta en x = 0. Ası que maxx∈[0,5] |f (2)(x)| = 1, y por lo tanto

|ET (f, 0,5)| ≤ 5 · (0,5)2

12· 1 = 0,1042

Para la regla de Simpson

|ES(f, h)| ≤ (b− a)h4

180maxx∈[0,5]

|f (4)(x)|

Similarmente, determinamos los puntos crıticos

f (5)(x) = (−15x+10x3−x5) exp(−x2/2) = 0 ⇒ x = 0 o x2 = 5±√

10.

En el segundo caso se tiene x ≈ 1,3556, x ≈ 2,8570. Comparando

f (4)(0) = 3

f (4)(1,3556) = −1,8549

f (4)(2,8570) = 0,3487

se encuentra que maxx∈[0,5] |f (4)(x)| = 3 y por lo tanto

|ES(f, h)| ≤ 5(0,5)4

180· 3 = 0,005208

(iv) Para la regla del trapecio M debe ser tal que

5(5/M)2

12≤ 5 · 10−9

de donde

M ≥ 5

√5

6104 = 45643,54

Entonces se necesita M = 45644. Para la regla del trapecio 2M debe ser talque

(5/2M)4

12≤ 5 · 10−9

de donde

M ≥ 5

2· 100

61/4= 159,73

Entonces se necesita M = 160.


Matlab

A continuacion las rutinas que implementan las reglas del trapecio y de Sim-pson compuestas.

function s = traprl (f, a, b, M)% Entrada - f es ei integrando introducido como cadena de caracteres ’f’% - a y b son los limites superior e inferior de integracion% - M es el numero de subintervalos% Salida - s es la suma de la regla trapezoidal (integral)h = (b - a) / M;s = 0;for k = 1:(M-1)

x = a + h * k;s = s + feval(f, x);

ends = h *(feval(f, a) + feval(f, b)) / 2 + h * s;function s = simprl (f, a, b, M)% Entrada - f es el integrando introducido como cadena de caracteres ’f’% - a y b son los limites superior e inferior de integracion% - M es el numero de subintervalos% Salida - s es la suma de la regla de simpsonh = (b - a) / (2 * M);s1 = 0;s2 = 0;for k = 1:M

x = a + h * (2 * k - 1);s1 = s1 + feval(f, x);

endfor k = 1:(M-1)

x = a + h * 2 * k;s2 = s2 + feval(f, x);

ends = h * (feval(f, a) + feval(f, b) + 4 * s1 + 2 * s2) / 3;

IV.5. Reglas Recursivas y Metodo de Romberg

IV.5.1. Reglas Recursivas

Trapecio. Para k = 0, 1, 2, . . ., el intervalo [a, b] se subdivide en Mk = 2k subin-tervalos de longitud hk = (b − a)/2k, y T (k) es el valor de la regla del trapeciocompuesta para estos 2k subintervalos. Las abscisas para T (k) son

xk,j = a+ jhk para j = 0, 1, . . . , 2k,

y los valores de la funcion son

fk,j = f(xk,j).

Los valores de las reglas del trapecio T (k) se pueden evaluar iterativamente de lasiguiente manera:

T0 =h

2(f(a) + f(b))

IV.5. REGLAS RECURSIVAS Y METODO DE ROMBERG IV.15

y para k ≥ 1

Tk =T (k − 1)

2+ hk

Mk−1∑j=1

fk,2j−1.

Prueba. Se verifica teniendo en cuenta que Mk = 2Mk−1 y que los nodos de Tk−1

son los nodos pares de Tk: es decir xk−1,j = xk,2j, y por lo tanto fk−1,j = fk,2j. Conesto, podemos escribir

Tk =hk

2

(fk,0 + 2

Mk−1∑j=1

fk,j + fk,Mk

)

=hk

2

fk,0 + 2

Mk2−1∑

j=1

fk,2j + fk,Mk

+ hk

Mk2∑

j=1

fk,2j−1

=1

2

hk−1

2

(fk−1,0 + 2

Mk−1−1∑j=1

fk−1,j + fk−1,Mk−1

)+ hk

Mk−1∑j=1

fk,2j−1

=Tk−1

2+ hk

Mk−1∑j=1

fk,2j−1.

Simpson. [a, b] se subdivide en 2k = 2M , k ≥ 1, subintervalos. Sk es el valorde la regla de Simpson compuesta con 2k subintervalos y hk = (b− a)/2k. Con uncalculo como el anterior, se obtiene la relacion

Sk =4Tk − Tk−1

3para k = 1, 2, 3, . . .

Es decir, por medio de esta formula se obtiene la cuadratura de Simpson con errorO(h4) a partir de la del trapecio con error O(h2).

Boole. [a, b] se subdivide en 2k = 4M , k ≥ 2, subintervalos. Bk es el valor de laregla de Boole compuesta con 2k subintervalos y hk = (b − a)/2k. De nuevo, conun calculo como el anterior, se obtiene la relacion

Bk =16Sk − Sk−1

15para k = 2, 3, . . .

Ası, se obtiene la cuadratura de Boole con error O(h6) a partir de la de Simpsoncon error O(h4).

Lo anterior es un ejemplo del metodo de Romberg el cual describimos a con-tinuacion.


IV.5.2. Metodo de Romberg

El metodo de Romberg se origina en la siguiente formula del error para la regadel trapecio:

ET (f, h) = a1h2 + a2h

4 + a3h6 + · · ·

Ya habıamos visto que este error es O(h2); lo crucial aquı es que cuando se de-rivan los terminos de mas alto orden no aparecen potencias impares, y que loscoeficientes ai son independientes de h. Ası que∫ b

a

f(x)dx = T (f, h) + a1h2 + a2h

4 + a3h6 + · · ·

Escribiendo la misma relacion para 2h se obtiene∫ b

a

f(x)dx = T (f, 2h) + 4a1h2 + 16a2h

4 + 64a3h6 + · · ·

Multiplicando por 4 la primera ecuacion y substrayendo la segunda se obtiene

(4− 1)

∫ b

a

f(x)dx = 4T (f, h)− T (f, 2h) + (4− 16)a2h4 + (4− 64)a3h

6 + · · ·

y ∫ b

a

f(x)dx =4T (f, h)− T (f, 2h)

4− 1+ b2h

4 + b3h6

con b2 = (4− 16)/(4− 1), b3 = (4− 64)/(4− 1), etc. De esta forma se ha obtenidouna nueva cuadratura con error O(h4). De acuerdo con las formulas recursivasderivadas antes, esta nueva regla es simplemente la regla de Simpson, la cualsabemos tiene error O(h4). Se puede continuar de la misma manera obteniendoreglas de error cada vez menor (la siguiente es la regla de Boole).

Generalizando lo anterior, se obtiene el siguiente procedimiento iterativo paracalcular una secuencia de cuadraturas Rk,j, j ≥ 1 y k ≥ j: primero

Rk,1 = Tk−1 para k ≥ 1

y para j ≥ 2 y k ≥ j se tiene la relacion

Rk,j =4j−1Rk,j−1 −Rk−1,j−1

4j−1 − 1

Esto corresponde a una tabla Rk,j calculada con las dependencias de la si-guiente manera:

IV.5. REGLAS RECURSIVAS Y METODO DE ROMBERG IV.17

Trapecio Simpson Boole · · ·

O(h2) O(h4) O(h6) O(h8) O(h10) O(h12) · · ·

R1,1

↖R2,1 ← R2,2

↖ ↖R3,1 ← R3,2 ← R3,3

↖ ↖ ↖R4,1 ← R4,2 ← R4,3 ← R4,4

↖ ↖ ↖ ↖R5,1 ← R5,2 ← R5,3 ← R5,4 ← R5,5

↖ ↖ ↖ ↖ ↖R6,1 ← R6,2 ← R6,3 ← R6,4 ← R6,5 ← R6,6

...

Se tiene que

Rk,1 = Tk−1 para k ≥ 1

Rk,2 = Sk−1 para k ≥ 2

Rk,3 = Bk−1 para k ≥ 3

Prueba. Veamos como se deduce la regla recursiva para Rk,j en general. Supon-

gamos que ya tenemos las reglas Rk,j−1 y Rk−1,j−1 las cuales tienen error O(h2(j−1)k )

y O(h2(j−1)k−1 ) respectivamente:∫ b

a

f(x)dx = Rk,j−1 + aj−1h2(j−1)k + ajh

2jk + · · ·∫ b

a

f(x)dx = Rk−1,j−1 + aj−1h2(j−1)k−1 + ajh

2jk−1 + · · ·

Usado hk−1 = 2hk, podemos reescribir la segunda ecuacion como∫ b

a

f(x)dx = Rk−1,j−1 + 22(j−1)aj−1h2(j−1)k + 22jajh

2jk · · ·

donde los puntos denotan terminos de mas alto orden. Multiplicando la primeraigualdad por 22(j−1) = 4j−1 y substrayendo la segunda se tiene que

(4j−1 − 1)

∫ b

a

f(x)dx =(4j−1Rk,j−1 −Rk−1,j−1

)+ bjh

2jk + · · ·

donde bj = −3 · 4j−1aj. Ası que∫ b

a

f(x)dx =4j−1Rk,j−1 −Rk−1,j−1

4j−1 − 1+ bjh

2jk + · · ·


y por lo tanto podemos definir la cuadratura

Rk,j =4j−1Rk,j−1 −Rk−1,j−1

4j−1 − 1

la cual tiene error O(h2jk ).

Ejemplo. Consideramos de nuevo la funcion f(x) = exp(−x2/2) en [0, 5].

>> [R Q err h]=romber(@(x) exp(-x.^2/2),0,5,4,5e-5)R =

2.50000932 0 0 0 01.35984699 0.97979288 0 0 01.25331998 1.21781098 1.23367885 0 01.25331289 1.25331052 1.25567716 1.25602633 01.25331327 1.25331340 1.25331359 1.25327608 1.25326529

Q =1.25326529

err =0.00276105

h =0.31250000

Matlab

function [R, quad, err, h] = romber (f, a, b, n, tol)% Entrada - f es el integrando introducido como cadena de caracteres ’f’% - a y b son los limites superior e inferior de integracion% - n es el numero maximo de filas en la tabla% - tol es la tolerancia% Salida - R es la tabla de Romberg% - quad es el valor de la cuadratura (integral)% - err es el error estimado% - h es el tama~no de paso mas peque~no utilizadoM = 1;h = b - a;err = 1;J = 0;R = zeros(4, 4);R(1, 1) = h * (feval(f, a) + feval(f, b)) / 2;while ((err > tol) & (J < n)) | (J < 4)

J = J + 1;h = h / 2;s = 0;for p = 1:M

x = a + h * (2 * p - 1);s = s + feval(f, x);

endR(J+1, 1) = R(J, 1) / 2 + h * s;M = 2 * M;for K = 1:J

R(J+1, K+1) = R(J+1, K) + (R(J+1, K) - R(J, K)) / (4 ^ K - 1);enderr = abs(R(J, J) - R(J+1, K+1));

endquad = R(J+1, J+1);

IV.6. CUADRATURA DE GAUSS-LEGENDRE IV.19

IV.6. Cuadratura de Gauss-Legendre

Las reglas que hemos estudiado hasta ahora tienen x0 = a y xn = b –por lo cualse dice que son cerradas– y los nodos estan distribuidos uniformemente. Existenreglas en las cuales esto no es cierto. En la cuadratura de Gauss-Legendre, con nnodos se logra grado de precision 2n− 1. Los nodos que se usan son las raıces delpolinomio de Legendre de grado n.

Regla de Gauss-Legendre para dos Puntos. La regla de Gauss-Legendrepara dos puntos esta dada por∫ 1

−1

f(x)dx ≈ G2(f) = f

(−1√

3

)+ f

(1√3

)con error

EG2(f) =f (4)(c)

135

para algun c ∈ [−1, 1], asumiendo f ∈ C4[−1, 1].

Esta formula se puede derivar a partir de la forma

G2(f) = w1f(α1) + w2f(α2)

tratando de obtener una cuadratura con el mayor grado de precision posible. Ası,procedemos con los calculos:∫ 1

−11dx = 2 G2(1) = w1 + w2∫ 1

−1xdx = 0 G2(x) = w1α1 + w2α2∫ 1

−1x2dx = 2

3G2(x

2) = w1α21 + w2α

22∫ 1

−1x3dx = 0 G2(x

3) = w1α31 + w2α

32

De donde se obtienen las ecuaciones

w1 + w2 = 2

w1α1 + w2α2 = 0

w1α21 + w2α

22 =

2

3w1α

31 + w2α

32 = 0

De la segunda ecuacion, se tiene w2α2 = −w1α1. Reemplazando esto en la ultimaresulta en

w1α1(α21 − α2

2) = 0

No puede ser que w1x1 = 0, por lo tanto α2 = −α1. Esto en la segunda ecuacionda w1 = w2, y por lo tanto la primera ecuacion implica que

w1 = w2 = 1.


Finalmente, la tercera ecuacion implica que

α1 = −α2 =1√3.

Los valores determinados corresponden a la expresion de G2 dada antes. El gradode precision es 3 porque∫ 1

−1x4dx = 2

5G2(x

4) =(− 1√

3

)4

+(

1√3

)4

= 29

Ejemplo. Usamos la cuadratura G2 para aproximar ln 5 =∫ 5

1dxx

. Primero debe-mos trasladar la integral al intervalo [−1, 1] con el cambio de variable x = 2u+ 3:∫ 5

1

dx

x=

∫ 1

−1

2du

2u+ 3.

Entonces

ln 5 =

∫ 1

−1

2du

2u+ 3≈ 2

2(−1/√

3) + 3+

2

2(1/√

3) + 3=

2 · 69− 4/3

=36

23= 1,565217 · · ·

mientras que el valor exacto es 1,609437912 · · ·

Regla de Gauss-Legendre para tres Puntos. La regla de Gauss-Legendrepara tres puntos esta dada por∫ 1

−1

f(x)dx ≈ G3(f) =1

9

(5f

(−√

3

5

)+ 8f(0) + 5f

(√3

5

))con error

EG3(f) =f (6)(c)

15750

para algun c ∈ [−1, 1], asumiendo f ∈ C6[−1, 1].

Ejemplo. Usamos la cuadratura G3 para aproximar ln 5 =∫ 5

1dxx

. Ya sabemos

ln 5 =

∫ 5

1

dx

x=

∫ 1

−1

2du

2u+ 3.

Entonces

ln 5 ≈ 1

9

(5

2

−2√

3/5 + 3+ 8

2

2 · 0 + 3+ 5

2

+2√

3/5 + 3

)=

2

5(5 · 0,68927172 + 8 · 0,33333333 + 5 · 0,21981919)

=2

9(3,44635859 + 2,66666667 + 1,09909596)

= 1,60269360 · · ·


IV.6.1. Polinomios de Legendre

Se puede usar el metodo anterior (coeficientes indeterminados) para obtenern nodos y coeficientes de tal manera que la cuadratura resultante tiene gradode precision 2n − 1. Alternativamente, esos mismos nodos resultan ser las raıcesde los llamados polinomios de Legendre. En el intervalo [−1, 1], el conjunto depolinomios forma un espacio vectorial con producto interior

〈P (x), Q(x)〉 =

∫ +1

−1

P (x)Q(x)dx

Se puede aplicar entonces el proceso de ortogonalizacion de Gram-Schmidt paraobtener una sucesion de polinomios ortogonales Pn(x), n = 0, 1, 2 . . .. Estos sonlos llamados polinomios de Legendre. Enunciamos el resultado que se obtiene.

Teorema IV.1 Los polinomios monicos de Legendre P0(x), P1(x), P2(x), . . . defi-nidos por la relacion de recurrencia

P0(x) = 1, P1(x) = x,

Pn+1(x) = xPn(x)− n2

4n2 − 1Pn−1(x)

satisfacen las siguientes propiedades: para cada n = 0, 1, 2, 3, . . .

(i) Pn(x) tiene n raıces reales diferentes en el intervalo [−1, 1]

(ii) cualquier polinomio P (x) de grado n se puede escribir como combinacionlineal de P0(x), P1(x), . . . , Pn(x), es decir, existen numeros reales c0, . . . , cn(unicos) tal que

P (x) =n∑

i=0

ciPi(x)

(iii) para todo m 6= n, Pn(x) y Pm(x) son ortogonales, es decir,∫ +1

−1

Pn(x)Pm(x)dx = 0

(iv) si Q(x) es de grado menor que n entonces∫ +1

−1

Pn(x)Q(x)dx = 0.

Como ejemplo, tenemos la siguiente tabla de los primeros polinomios y susraıces:


n Pn(x) raıces0 11 x 0

2 x2 − 13

±√

13

= ±√

33

3 x3 − 35x 0,±

√35

=√

155

4 x4 − 67x+ 30

35±√

17

(3± 2

√65

)= ±√

525±70√

30

35

Note que las raıces de los polinomios de grado 2 y 3 son las mismas obtenidasantes usando el metodo de coeficientes indeterminados.

IV.6.2. Metodo de Gauss-Legendre

El metodo de Gauss-Legendre con n nodos para la integral∫ +1

−1f(x)dx es la

cuadratura

Gn[f ] =n∑

i=1

wif(xi)

donde x1, x2, . . . , xn son las n raıces de Pn(x), y los pesos wi son los correspondien-tes a la integral exacta del polinomio de interpolacion para los datos (xi, f(xi)),i = 1, 2, . . . , n. Mas precisamente, recordemos que para estos datos la interpolacionde Lagrange (un polinomio de grado n− 1) esta dada por

I(x) =n∑

i=1

f(xi) · Ln−1,i(x),

donde

Ln,i =∏j 6=i

x− xj

xi − xj

.

Entonces ∫ +1

−1

I(x)dx =n∑

i=1

f(xi)

∫ +1

−1

Ln−1,i(x)dx,

y por lo tanto

wi =

∫ +1

−1

Ln−1,i(x)dx.

Esto da un metodo para calcular los wi de la cuadratura. Equivalentemente, loswi se pueden obtener con el metodo de coeficientes indeterminados para los nodosx1, . . . , xn y los polinomios 1, x, x2, . . . , xn−1, de tal manera que la cuadratura esexacta para los polinomios 1, x, x2, . . . , xn−1. Se obtienen las n ecuaciones, parak = 0, 1, 2, . . . , n− 1, con las n incognitas w1, w2, . . . , wn:

n∑i=1

wi · xki =

∫ +1

−1

xkdx =1− (−1)k+1

k + 1.


Este metodo se ha implementado en el programa Matlab que se incluye adelante.

Ahora justificamos el grado de precision del metodo de Gauss-Legendre.

Teorema IV.2 El metodo de Gauss-Legendre es exacto para todo polinomio degrado menor que 2n.

Prueba. Sea P (x) un polinomio de grado menor que 2n. Si el grado de P (x) esmenor que n entonces la cuadratura es exacta (por la construccion de los wi). SiP (x) tiene grado al menos n y menor que 2n, entonces tomando el cociente conPn(x) podemos escribir

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

donde los grados del cociente Q(x) y del residuo R(x) son menores que n. Por lotanto, por el caso ya considerado, la cuadratura para R es exacta, es decir∫ +1

−1

R(x)dx = Gn[R].

Por otra parte∫ +1

−1

P (x)dx =

∫ +1

−1

Q(x)Pn(x)dx+

∫ +1

−1

R(x)dx =

∫ +1

−1

R(x)dx

porque∫ +1

−1Q(x)Pn(x)dx = 0 por la propiedad (iv) de los polinomios de Legendre

en el teorema anterior. Ademas puesto que x1, x2, . . . , xn son las raıces de Pn(x),es decir Pn(xi) = 0, se tiene que

P (xi) = Q(x)Pn(xi) +R(xi) = R(xi),

y por lo tanto las cuadraturas de P (x) y de R(x) son iguales, Gn[P ] = Gn[R].Ası que ∫ +1

−1

P (x)dx =

∫ +1

−1

R(x)dx = Gn[R] = Gn[P ].

Hemos obtenido que la cuadratura es exacta para P (x) donde P (x) es cualquierpolinomio de grado menor que 2n.

Finalmente, en general la integral de interes es sobre un intervalos [a, b]. Parasu evaluacion, se usa un cambio de variable

x =a(1− t) + b(1 + t)

2Con esto, ∫ b

a

f(x)dx =b− a

2

∫ +1

−1

f

(a(1− t) + b(1 + t)

2

)dt

y por lo tanto la cuadratura de f en [a, b] con n nodos esta dada por

Gn[f ] =b− a

2·

n∑i=1

wi · f(a(1− xi) + b(1 + xi)

2

)donde los xi son las raıces en [−1, 1] como antes, y los coeficietes wi son los mismosque antes.


IV.6.3. Programa Matlab

El siguiente programa Matlab calcula las raıces del polinomio de Legendrenecesario y los coeficiente de la cuadratura. Por supuesto, si esta funcion se vaa usar repetidamente serıa mas eficiente precalcular estos, tenerlos disponibles ysimplemente usarlos. El calculo de las raıces se basa en la observacion de que elpolinomio de Legendre Pn(x) es igual al polinomio caracterıstico de la siguientematriz, llamada matriz de Jacobi:

0 1/√

3 0 0 0 0 · · · 0 0 01/√

3 0 2/√

15 0 0 0 · · · 0 0 00 2/

√15 0 3/

√35 0 0 · · · 0 0 0

0 0 3/√

35 0 4/√

65 0 · · · 0 0 0...

......

......

.... . .

......

...0 0 0 0 0 0 · · · n−2√

4(n−2)2−10 n−1√

4(n−1)2−1

0 0 0 0 0 0 · · · 0 n−1√4(n−1)2−1

0

Esto se puede comprobar de la ecuacion de recurrecia de los polinomios. La eva-luacıon de los coeficientes se basa en el metodo de coeficientes indeterminados.

function [quad, raices, nodos, coefs] = GaussLegendre( f, a, b, N)% Entrada - f es el integrando introducido como una cadena de caracters ’f’% - a y b limites superior e inferior de integracion% - N es el numero de nodos en la cuadratura% Salida - quad es el valor de cuadratura% - raices es el vector de raices del polinomio de Legendre N% - nodos es el vector de nodos de la cuadratura e [a,b]% - coefs es el vector de pesos de la cuadratura

% evaluacion de las N raices del polinomio de Legendre P_N de grado N,% usando el hecho de que P_N es el polinomio caracteristico de la% matriz de Jacobi J

v = feval(@(n) n./sqrt(4*n.^2-1), 1:N-1);J = diag(v,-1)+diag(v,+1);raices = eig(J);

% evaluacion de los nodos por medio de transformacion lineal

nodos = (a*(1-raices)+b*(1+raices))/2;

% evaluacion de los coeficientes de la cuadratura por el metodo de% coeficientes indeterminados

V = (vander(raices))’;z = (feval(@(k) (1-(-1).^k)./k, N:-1:1))’;coefs = V\z;

% evaluacion de la cuadratura

quad = ((b-a)/2) * sum(coefs.*feval(f, nodos));

end


Ejemplo. Usamos la rutina para evaluar∫ 5

0exp(−x2/2)dx con 8 nodos:

>> [quad, raices, nodos, coefs] = GaussLegendre(@(x) exp(-x.^2/2),0,5,8)

quad =1.253315160546044

raices =-0.960289856497536-0.796666477413627-0.525532409916329-0.1834346424956500.1834346424956500.5255324099163290.7966664774136270.960289856497536

nodos =0.0992753587561590.5083338064659331.1861689752091782.0414133937608752.9585866062391263.8138310247908224.4916661935340664.900724641243841

coefs =0.1012285362903750.2223810344533770.3137066458778840.3626837833783670.3626837833783580.3137066458778900.2223810344533730.101228536290377

Cap tulo IV Integraci on Num erica - Sede Medellinearamosn/cursos/... · 2013. 10. 28. · Cap tulo...

Documents

Transcript of Cap tulo IV Integraci on Num erica - Sede Medellinearamosn/cursos/... · 2013. 10. 28. · Cap tulo...