¡Descarga Integración numérica y más Apuntes en PDF de Cálculo solo en Docsity!
C´alculo Num´erico I
Integraci´on Num´erica
La integral resuelve el problema de calcular el ´area bajo la gr´afica de una funci´on positiva definida sobre un intervalo cerrado. El c´alculo elemental de funciones de una variable real proporciona un m´etodo elegante de calcular la integral de una funci´on. El teorema funda- mental del c´alculo nos dice que el problema de calcular la integral de una funci´on continua se reduce al de buscar una segunda funci´on cuya derivada sea la funci´on dada, es decir una primitiva de ella. Sin embargo el problema de hallar una primitiva de una funci´on dada puede resultar muy dif´ıcil si no imposible. De hecho sabemos que existen funciones elementales —es decir, combinaciones alge- braicas de funciones trigonom´etricas y logar´ıtmicas y sus inversas— cuyas primitivas no son expresables de esta forma (p. ej. e−x 2 ). Por esta raz´on es por la que estudiamos m´etodos num´ericos que apro- ximen el valor de la integral buscada. Ya la definici´on de integral de Riemann proporciona un m´etodo de aproximaci´on num´erica: las sumas de Riemann. Sin embargo su convergencia es muy lenta y no resultan ´utiles para obtener resultados pr´acticos. Los m´etodos num´ericos que vamos a estudiar consisten en sustituir la funci´on dada por una aproximaci´on suya y tomar como valor de la integral de la funci´on el valor de la integral de su aproximada.
Veremos en primer lugar los resultados que se obtienen aproxi- mando la funci´on por medio de un polinomio interpolador con es- pecial ´enfasis en los casos lineal (regla del trapecio) y cuadr´atico (regla de Simpson). A continuaci´on estudiaremos las cuadraturas de Gauss.
18 Regla del Trapecio
Si para calcular el valor aproximado de la funci´on f en el inter- valo [a, b ] sustituimos dicha funci´on por el polinomio lineal que la interpola con nodos en los extremos del intervalo obtenemos
T f =
∫ (^) b
a
(x − a)f (b) + (b − x)f (a) b − a
= (b − a)
f (b) + f (a) 2
¿Cu´al ser´a el error cometido?
E(T f ) = T f − I(f ) = −
∫ (^) b
a
(x − a)(x − b)f [a, b, x]dx
Por el teorema del valor medio generalizado del c´alculo integral
∃ξ ∈ (a, b) : E(T f ) = −f [a, b, ξ]
∫ (^) b
a
(x − a)(x − b)dx
y por tanto, si f ∈ C^2 ([a, b]),
∃ζ ∈ (a, b) : E(T f ) =
(b − a)^3 12
f ′′(ζ).
Si b − a no es peque˜no la regla del trapecio no ser´a muy ´util para calcular I(f ). En ese caso podr´ıamos aplicarla dividiendo antes el intervalo [a, b ] en un cierto n´umero n de subintervalos de lon- gitud h = (b − a)/n con extremos xj = a + jh, j = 0, 1 ,... , n. Tendr´ıamos entonces
I(f ) =
∑^ n
j=
∫ (^) xj
xj− 1
f (x)dx =
∑^ n
j=
h 2
[f (xj− 1 ) + f (xj )] −
h^3 12
f ′′(ξj )
= h[^12 f (x 0 ) + f (x 1 ) + · · · + f (xn− 1 ) + 12 f (xn)] −
h^3 12
∑^ n
j=
f ′′(ξj ).
Llamamos regla del trapecio compuesta a regla del trapecio compuesta
Tnf =
b − a
n
[^12 f (x 0 ) + f (x 1 ) + · · · + f (xn− 1 ) + 12 f (xn)].
El error que se comete al tomar Tnf en vez de If es
E(Tnf ) = Tn(f ) − I(f ) =
h^3 12
∑^ n
j=
f ′′(ξj )
h^2 12
(b − a)
n
∑^ n
j=
f ′′(ξj ) =
h^2 12
(b − a)f ′′(ξ), ξ ∈ (a, b).
Este error puede estimarse asint´oticamente de la siguiente manera:
lim n→∞
E(Tnf ) h^2
= lim n→∞
h 12
∑^ n
j=
f ′′(ξj )
cantidad esta ´ultima que se puede interpretar como una suma de Riemann de forma que
lim n→∞
E(Tnf ) h^2
∫ (^) b
a
f ′′(x)dx =
f ′(b) − f ′(a) 12
Para determinar el error en S(f ) integraremos el error en el poli- nomio interpolador de segundo grado
E(Sf ) = −
∫ (^) b
a
(x − a)(x − c)(x − b)f [a, c, b, x]dx.
Para calcular esta integral no podemos aplicar directamente el teo- rema del valor intermedio del c´alculo integral ya que ahora la funci´on (x − a)(x − c)(x − b) cambia de signo en el punto c. Sorteamos este problema por medio de la funci´on
w(x) =
∫ (^) x
a
(t − a)(t − c)(t − b)dt,
primitiva de la anterior, que verifica las propiedades: w(a) = 0 = w(b) y ∀x ∈ (a, b), w(x) > 0. Calculamos entonces la integral en E(Sf ) integrando por partes:
∫ (^) b
a
(P 2 (x) − f (x))dx = −
∫ (^) b
a
w′(x)f [a, b, c, x]dx
= −(w(x)f [a, b, c, x])
∣ba +
∫ (^) b
a
w(x)
d dx
f [a, b, c, x]dx
∫ (^) b
a
w(x)
d dx
f [a, b, c, x]dx.
Y a esta ´ultima integral s´ı que puede aplicarse el teorema del valor medio: (^) dxd f [a, b, c, x] = f [a, b, c, x, x]
= f [a, b, c, ξ, ξ]
∫ (^) b
a
w(x)dx.
Dado que, si escribimos, h = (b − a)/2,
∫ (^) b
a
w(x)dx =
∫ (^) b
a
∫ (^) x
a
(t − a)(t − b)(t − c)dtdx
∫ (^) b
a
∫ (^) x−a−h
−h
(u + h)(u − h)u dudx
∫ (^) h
−h
∫ (^) y
−h
(u^3 − h^2 u) dudy
= 4 h^5 / 15.
Obtenemos
E =
f (4)(ξ) 24
b − a 2
)^5.
Nuevamente, si h = (b − a)/2 es grande el error ser´a inadmisible. Como hicimos en el caso de la regla del trapecio podemos aplicar la regla de Simpson a trozos sobre el intervalo [a, b ]:
sea h = (b − a)/ 2 n, xk = a + kh, jk = 0, 1 ,... , 2 n,
∫ (^) b
a
f =
∑^ n
j=
∫ (^) x 2 j
x 2 j− 2
f
∑^ n
j=
h 3
[f (x 2 j− 2 ) + 4f (x 2 j− 1 ) + f (x 2 j )] −
∑^ n
j=
h^5 90
f IV^ (ξj ).
Llamamos regla de Simpson compuesta a regla de Simpson compuesta
Snf =
b − a
6 n
[f (x 0 ) + 4f (x 1 ) + 2f (x 2 ) + 4f (x 3 )+
· · · + 2f (x 2 n− 2 ) + 4f (x 2 n− 1 ) + f (x 2 n)]
El error que se comete al tomar Snf en vez de If es
E(Snf ) =
h^5 90
∑^ n
j=
f iv(ξj ) =
h^4 90
(b−a)
2 n
∑^ n
j=
f iv(ξj ) =
h^4 180
(b−a)f iv(ξ).
Tambi´en ahora podemos estimar el error por medio de
∑^ n
j=
h^5 90
f iv(ξj ) =
h^4 180
∑^ n
j=
2 hf iv(ξj ) ≈
h^4 180
∫ (^) b
a
f iv(ξ)dξ =
h^4 180
(f ′′′(b)−f ′′′(a)).
Si se corrige la regla de Simpson compuesta por medio de este error asint´otico se obtiene la regla de Simpson corregida: regla de Simpson corregida
CSn(f ) = Snf −
(b − a)^4
2880 n^4
[f ′′′(b) − f ′′′(a)].
Ejemplo 19.1. Tratando ahora el ejemplo anterior
∫ (^1)
0
4 dx 1 + x^2
con la regla de Simpson compuesta, y observando que 96 es una cota superior para | (^) dxd1+^4 x 2 | en el intervalo [0, 1], se obtiene que
|E(Snf )| ≤
90 n^4
buscamos n tal que 96 90 n^4
≤ 5 · 10 −^3
- Si n es par y f ∈ Cn+2[a, b], ∃ξ ∈ (a, b) tal que
Inf − If = Cnhn+3f (n+2)(ξ)
donde
Cn = −
(n + 2)!
∫ (^) n
0
t^2 (t − 1)... (t − n)dt.
- Si n es impar y f ∈ Cn+1[a, b], ∃ξ ∈ [a, b] tal que
Inf −
∫ (^) b
a
f = Cnhn+2f (n+1)(ξ)
donde
Cn = −
(n + 1)!
∫ (^) n
0
t(t − 1)... (t − n)dt.
21 Reglas Gaussianas
Las reglas de Newton–Cotes garantizan la integraci´on exacta de polinomios hasta un determinado grado: el error en las reglas im- pares (trapecio, tres octavos,...) depende de la derivada de un orden superior (segunda, cuarta,...) y por tanto la regla integra exacta- mente polinomios hasta el mismo grado que indica la regla (primero, tercero,...); el error en las reglas pares (punto medio, Simpson, Boole,...) depende de la derivada de orden dos unidades por encima de la regla (segunda, cuarta, sexta,...) y por tanto integra exacta- mente polinomios polinomios hasta un grado por encima del orden de la regla (primero, tercero, quinto,...).
A la vista de este comportamiento vamos a buscar reglas de in- tegraci´on num´erica de la forma
wif (xi) que integren todos los polinomios hasta un determinado grado, el mayor posible para el n´umero de sumandos elegidos para la regla. Obtendremos as´ı las llamadas reglas Gaussianas o cuadraturas de Gauss. reglas Gaussianas
Si utilizamos un s´olo nodo la regla resultante integrar´a exactamente cuadraturas de Gauss polinomios de grado menor o igual a uno: (normalizamos el c´alculo al intervalo [− 1 , 1]). Sea x 0 el nodo y w 0 el peso correspondiente, peso
∫ (^1)
− 1
dx = w 0 ,
− 1
xdx = w 0 x 0.
De aqu´ı resulta x 0 = 0, w 0 = 2.
Si utilizamos dos nodos podemos integrar exactamente todos los polinomios de grado menor o igual que tres:
I 2 f = w 1 f (x 1 ) + w 2 f (x 2 )
I 2 (1) = w 1 + w 2 = 2 I 2 (x) = w 1 x 1 + w 2 x 2 = 0 I 2 (x^2 ) = w 1 x^21 + w 2 x^22 = (^23) I 2 (x^3 ) = w 1 x^31 + w 2 x^32 = 0
Este sistema no lineal resulta tener soluci´on ´unica. Esta verifica´ w 1 = w 2 = 1 y x 1 = −x 2 = −
3 /3. As´ı
I 2 (f ) = f (−
) + f (
En general y por el mismo procedimiento, utilizando n nodos pode- mos integrar exactamente todos los polinomios de grado ≤ 2 n − 1. Para obtener los nodos y los pesos de In debemos resolver el sistema no lineal 2 = w 1 + · · · + wn 0 = w 1 x 1 + · · · + wnxn 2 3 =^ w^1 x
2 1 +^ · · ·^ +^ wnx 2 n
............ 2 2 n− 1 =^ w^1 x
2 n− 2 1 +^ · · ·^ +^ wnx 2 n− 2 n 0 = w 1 x^21 n −^1 + · · · + wnx^2 nn−^1
La resoluci´on de este sistema no lineal no es en general sencilla, de hecho no est´a claro a priori que el sistema tenga soluci´on o que, en caso de existir, ´esta sea ´unica.
Abordamos el problema de su resoluci´on por otros medios. Su- pongamos primeramente que existen los xi y que son distintos y est´an en el intervalo [0, 1]. El polinomio Ln(x) = (x − x 1 )(x − x 2 ) · · · (x − xn) debe tener integral exacta por medio de la regla de Gauss y por tanto esta debe ser cero, ya que Ln(xi) = 0. Si Ln(x) = xn^ + a 1 xn−^1 + · · · + an− 1 x + an podemos determinar los coeficientes ai por medio del sistema lineal de ecuaciones
− 1
xkLn(x)dx k = 1,... , n.
De hecho se obtienen dos sistemas lineales separados, uno que con- tiene solamente a los coeficientes de ´ındice par, y el otro solamente a los de ´ındice impar. Uno de ellos ser´a homog´eneo con soluci´on trivial, por tanto el polinomio obtenido tendr´a solamente t´erminos cuyo grado tenga la misma paridad que n. Los ceros de este poli- nomio son los nodos xi buscados. En general se obtendr´an por alg´un m´etodo de resoluci´on de ecuaciones no lineales (Newton,...).
Los polinomios Ln reciben el nombre de polinomios m´onicos de Legendre y han sido estudiados extensivamente. polinomios m´onicos de Legendre
