Unha introdución á integración numérica: fórmulas de Gauss

Traballo Fin de Grao UNHA INTRODUCIÓN Á INTEGRACIÓN NUMÉRICA. FÓRMULAS DE GAUSS. Antía Álvarez González Curso 2021-2022 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA GRAO DE MATEMÁTICAS Traballo Fin de Grao UNHA INTRODUCIÓN Á INTEGRACIÓN NUMÉRICA. FÓRMULAS DE GAUSS. Antía Álvarez González Xullo, 2022 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA Traballo proposto Área de Coñecemento: Matemática Aplicada Título: Unha introdución á integración numérica. Fórmulas de Gauss. Breve descrición do contido Realízase unha introdución aos métodos de cuadratura de tipo interpolatorio-polinómico, para despois profundar na análise destas fórmulas ata chegar ás denominadas fórmulas de Gauss. iii Índice Resumo viii Introdución xi 1. Preliminares 1 1.1. Fórmulas de tipo interpolatorio-polinómico . . . . . . . . . . . . . . . . . . . . . 1 1.2. Fórmulas de Newton-Cotes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.2.1. Expresión do erro nas fórmulas de Newton-Cotes . . . . . . . . . . . . . . 8 1.2.2. Obtención das fórmulas de Newton-Cotes . . . . . . . . . . . . . . . . . . 9 1.3. Fórmulas de Newton-Cotes compostas . . . . . . . . . . . . . . . . . . . . . . . . 12 2. Fórmulas de Gauss 15 2.1. Polinomiosortogonais ................................. 16 2.2. Fórmulas de cuadratura de Gauss . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.2.1. Máximo grao de exactitude . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.2.2. Unicidade das fórmulas de Gauss . . . . . . . . . . . . . . . . . . . . . . . 26 2.2.3. Estimación do erro nas fórmulas de Gauss . . . . . . . . . . . . . . . . . . 27 2.3. Cálculo das fórmulas de cuadratura gaussianas . . . . . . . . . . . . . . . . . . . 27 3. Integración de Romberg 33 3.1. Extrapolación de Richardson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 v vi ÍNDICE 3.2. IntegracióndeRomberg ................................ 34 Bibliografía 39 2 1. Preliminares modo, reemplázase o cálculo da integral I(f) = Zb a f(x)dx polo de Zb a pn(x)dx , de modo que se considera a seguinte aproximación para a integral: Zb a f(x)dx ≈Zb a pn(x)dx. A principal vantaxe deste tipo de fórmulas é a que a integración de polinomios é un proceso que, polo xeral, non require de cálculos demasiado complexos. Nótese ademais que o erro de cuadratura, que xa se deniu na expresión (1.2), se pode reescribir nas fórmulas de tipo interpolatorio-polinómico como segue: En(f) = I(f)−In(f) = Zb a (f(x)−pn(x))dx, é dicir, como a integral do erro cometido na interpolación da función f . Fixados n+ 1 nodos calesquera situados sobre o intervalo [a, b] , tomarase como polinomio de interpolación pn o polinomio interpolatorio de Lagrange relativo a eses n+ 1 nodos. Denición 1.2. O polinomio interpolatorio de Lagrange relativo a n+ 1 puntos , calcúlase mediante pn(x) = n X i=0 f(xi)li(x), onde li(x) é o polinomio característico de Lagrange asociado ao nodo xi , que se dene como li(x) = n Y k=0 k6=i x−xk xi−xk . (1.3) Para cada función f e conxunto de n+ 1 puntos, existe un único polinomio interpolatorio de Lagrange que cumpre esta denición. Pódese ver a proba deste resultado, por exemplo, en Stoer e Bulirsch (2010). A partir deste polinomio, a fórmula de cuadratura para o cálculo da integral I(f) pódese reescribir como In(f) = Zb a n X i=0 f(xi)li(x)dx = n X i=0 f(xi)Zb a li(x)dx. (1.4) Nótese que esta fórmula ten a forma In(f) = n X i=0 αif(xi) , onde xi son os nodos de cuadratura e αi=Zb a li(x)dx son os pesos para cada i= 0, ... , n , polo que claramente é do tipo intepolatoriopolinómico. De acordo con esto, obtense un método para o cálculo dos pesos de cuadratura a través da integración dos polinomios característicos de Lagrange asociados a cada un dos nodos. Móstrase a continuación un exemplo que ilustra proceso. 1.1. Fórmulas de tipo interpolatorio-polinómico 3 Exemplo 1.3. Cálculo dos pesos de cuadratura empregando os polinomios característicos de Lagrange para unha fórmula de cuadratura de dous nodos nun intervalo [a, b] . Considéranse os nodos de cuadratura x0=a e x1=b . Polo tanto, os pesos calcúlanse como α0=Zb a l0(x)dx =Zb a x−x1 x0−x1 dx =Zb a x−b a−bdx = =1 a−bx2 2−bxb a =b−a 2 α1=Zb a l1(x)dx =Zb a x−x0 x1−x0 dx =Zb a x−a b−adx = =1 b−ax2 2−axb a =b−a 2 de modo que α0=α1=b−a 2 .  A seguinte denición será de gran importancia para determinar como de exacta é unha aproximación, e poder comparar así as distintas clases de fórmulas de tipo interpolatorio-polinómico que se verán posteriormente. Denición 1.4. O grao de exactitude ou de precisión dunha fórmula de cuadratura defínese como o maior número enteiro positivo r para o cal a fórmula de cuadratura é exacta para todo polinomio de grao ≤r , é dicir, para o cal En(xk) = 0 para cada k= 0,1, ... , r . Teorema 1.5. Unha fórmula de cuadratura da forma n X i=0 αif(xi) de n+ 1 puntos é de tipo interpolatorio-polinómico se e só se ten, polo menos, grao de exactitude n . Demostración. ( ⇒ ) Dados n+1 puntos situados no intervalo [a, b] e unha fórmula de cuadratura In(f) = n X i=0 αif(xi) de tipo interpolatorio-polinómico, xa se viu que os seus pesos se poden calcular mediante αi=Zb a li(x)dx, i = 0, ... , n, sendo li o polinomio característico de Lagrange asociado ao nodo xi para cada i= 0, ... , n . Denotarase Pn como o conxunto dos polinomios de grao menor ou igual que n . Como todo polinomio p∈Pn de grao n coincide co seu polinomio de interpolación relativo aos puntos x0, ... , xn , tense que p(x) = n X i=0 p(xi)li(x), 4 1. Preliminares e polo tanto integrando esto, como αi=Zb a li(x)dx , queda Zb a p(x)dx = n X i=0 αip(xi), polo que a fórmula é exacta para o polinomio p∈Pn . ( ⇐ ) Para probar o recíproco, suporase que a fórmula é exacta para todo polinomio de Pn . En particular, ten que selo para lk∈Pn de grao n , e así: Zb a lk(x)dx = n X i=0 αilk(xi) = n X i=0 αiδki =αk, k = 0, ... , n, logo a fórmula é de tipo interpolatorio-polinómico, xa que os pesos se calculan como se viu en (1.4). Nótese que a través deste resultado, se proporciona un método que permite a obtención dos pesos de cuadratura sen calcular a integral dos polinomios característicos de Lagrange, Zb a li(x)dx , para cada i= 0, ... , n . Dada unha fórmula de cuadratura de n+ 1 nodos, polo teorema previo, esta terá grao de exactitude n , e polo tanto será exacta para todo monomio xk tal que 0≤k≤n . Polo tanto, se x0, ... , xn son os n+ 1 nodos, debe satisfacerse o seguinte sistema linear de n+ 1 ecuacións                      α0+α1+... +αn=Zb a 1dx α0x0+α1x1+... +αnxn=Zb a xdx . . . α0xn 0+α1xn 1+... +αnxn n=Zb a xndx (1.5) que tamén se pode expresar matricialmente como V α =d, (1.6) onde α é o vector de pesos, d= (d0, d1, ... , dn)T, con dk=Zb a xkdx, k = 0, ... , n, e V é a matriz de Vandermonde trasposta: V=          1 1 ··· 1 x0x1··· xn x2 0x2 1··· x2 n . . .. . .. . . xn 0xn 1··· xn n          . 1.2. Fórmulas de Newton-Cotes 5 O sistema (1.6) terá solución única se o determinante da matriz V é non nulo, é dicir, se se cumpre det(V) = Y 1≤j<i≤n (xi−xj)6= 0, o cal debe suceder porque os nodos de cuadratura escóllense todos distintos dous a dous. Empregando os mesmos parámetros ca no Exemplo 1.3, móstrase a cotinuación o proceso do cálculo dos pesos de cuadratura empregando o novo método proporcionado polo sistema (1.5). Exemplo 1.6. Cálculo dos pesos de cuadratura a través do sistema (1.5) . Os nodos de cuadratura son x0=a e x1=b . O sistema linear mencionado queda neste caso como        α0+α1=Zb a 1dx =b−a α0a+α1b=Zb a x dx =b2−a2 2 Resolvéndoo, vese que os pesos obtidos son α0=α1=b−a 2. que coinciden cos do Exemplo 1.3.  1.2. Fórmulas de Newton-Cotes Dentro das fórmulas de cuadratura de tipo interpolatorio-polinómico, estudarase primeiro a familia das fórmulas de Newton-Cotes. Caracterízase porque se empregan n+1 nodos distribuídos de xeito equidistante sobre o intervalo [a, b] . Dado un paso h determinado (que se dene como a distancia entre dous nodos), calcúlanse como xi=x0+ih , para cada i= 0, ... , n . Existen dous tipos de fórmulas de Newton-Cotes dependendo de como se escollan os extremos do intervalo: Pechadas: inclúense os extremos do intervalo de integración, polo que se toma x0=a, xn= b e h= (b−a)/n para n≥1 . Abertas: non se inclúen os extremos do intervalo de integración, polo que x0=a+h, xn= b−h e h= (b−a)/(n+ 2) para n≥0 . Proposición 1.7. Nas fórmulas de Newton-Cotes, os pesos de cuadratura só dependen dos valores de h e de n . Demostración. Introdúcese o cambio de variable x= Ψ(t) = x0+th . No caso pechado considéranse os nodos xk=a+kh para cada k= 0, ... , n , e polo tanto x−xk xi−xk =a+th −(a+kh) a+ih −(a+kh)=t−k i−k. 6 1. Preliminares Entón, como para n≥1 li(x) = n Y k=0 k6=i t−k i−k=ϕi(t),0≤i≤n, os pesos de cuadratura pódense expresar como αi=Zb a li(x)dx =Zn 0 ϕi(t)hdt =hZn 0 ϕi(t)dt. (1.7) Polo tanto, a fórmula de cuadratura queda como In(f) = h n X i=0 wif(xi), wi=Zn 0 ϕi(t)dt. (1.8) Para o caso aberto, considéranse os nodos xk=a+h(k+ 1) , e denotando os extremos do intervalo [a, b] por x−1=a e xn+1 =b , e co mesmo cambio de variable que antes, os pesos son αi=hZn+1 −1 ϕi(t)dt . Así, a fórmula de cuadratura neste caso quedará como In(f) = h n X i=0 wif(xi), wi=Zn+1 −1 ϕi(t)dt. (1.9) Tendo en conta que ϕi(t) = n Y k=0 k6=i t−k i−k dedúcese, por como están denidos os coecientes wi en (1.8) e (1.9) para cada i= 0, ... , n , que no cálculo de In(f) únicamente se teñen en conta os valores de h e de n . Obsérvese que a proba deste resultado proporciona un novo método para o cálculo dos pesos de cuadratura das fórmulas de Newton-Cotes, a través da expresión (1.7). O exemplo que se dá a continuación ilustra este proceso. Débese destacar ademais o fácil que resulta o cálculo dos pesos mediante este método en comparación cos vistos anteriormente, polo que tamén se compara co método que deriva do sistema linear (1.5). Exemplo 1.8. Cálculo dos pesos de cuadratura dunha fórmula de Newton-Cotes pechada con n=2. (a) Mediante a forma proporcionada polo sistema (1.5) . Neste caso h=b−a 2 , polo que os nodos de cuadratura son x0=a, x1=a+b 2 e x2=b . Considérase polo tanto o seguinte sistema de ecuacións                α0+α1+α2=Zb a 1dx =b−a α0a+α1 a+b 2+α2b=Zb a x dx =b2−a2 2 α0a2+α1 (a+b)2 4+α2b2=Zb a x2dx =b3−a3 3 1.2. Fórmulas de Newton-Cotes 7 con solución α0=α2=b−a 6 α1=2(b−a) 3. (b) Forma proporcionada pola expresión (1.7) . Recórdese que n= 2 e que h=b−a 2 , polo que os pesos se calculan como α0=hZ2 0 ϕ0(t)dt =hZ2 0 t−1 0−1 t−2 0−2dt =h 2Z2 0 (t−1)(t−2)dt = =h 2t3 3−3t2 2+ 2t2 0 =b−a 4 2 3=b−a 6 α1=hZ2 0 ϕ1(t)dt =hZ2 0 t−0 1−0 t−2 1−2dt =−hZ2 0 t(t−2)dt = =−ht3 3−t22 0 =−b−a 2−4 3=2(b−a) 3 α2=hZ2 0 ϕ2(t)dt =hZ2 0 t−0 2−1 t−1 2−1dt =h 2Z2 0 t(t−1)dt = =h 2t3 3−t2 22 0 =b−a 4 2 3=b−a 6 Claramente, os pesos obtidos coinciden cos de (a) .  Volvendo á Proposición 1.7, o feito de que os pesos de cuadratura das fórmulas de Newton- Cotes só dependan de h e de n , implica que, tanto a función f como o intervalo de integración que se consideren, non teñan relevancia na obtención destes coecientes. Isto implica que ditos pesos poidan ser calculados previamente á aplicación da fórmula de cuadratura. De feito, os pesos xa están predeterminados para moitos valores de n . Móstranse nas seguinte táboa os pesos αi nas fórmulas de Newton-Cotes para diferentes valores de n , primeiro no caso pechado e posteriormente no aberto. n α0α1α2α3α4 1 b−a 2 b−a 2 2 b−a 6 2(b−a) 3 b−a 6 3 b−a 8 3(b−a) 8 3(b−a) 8 b−a 8 4 7(b−a) 90 16(b−a) 45 6(b−a) 45 16(b−a) 45 7(b−a) 90 Cadro 1.1: Pesos nas fórmulas de Newton-Cotes pechadas 8 1. Preliminares n α0α1α2α3 0 b−a 1 b−a 2 b−a 2 2 2(b−a) 3−b−a 3 2(b−a) 3 3 11(b−a) 24 b−a 24 b−a 24 11(b−a) 24 Cadro 1.2: Pesos nas fórmulas de Newton-Cotes abertas Nótese que os pesos son simétricos, é dicir, que, no caso pechado tense que αi=αn−i para i= 0, ... , n −1 e no caso aberto, αi=αn−i para i= 0, ... , n . Ademais, tamén existe a presenza de pesos negativos. No caso aberto, a partir de n= 2 e no pechado, a pesar de que non se mostra no Cadro 1.1, a partir de n= 8 . Isto repercute de forma negativa nas fórmulas de cuadratura, xa que dá lugar a inestabilidades. 1.2.1. Expresión do erro nas fórmulas de Newton-Cotes Xa se viu no Teorema 1.5 que ao aproximar unha integral I(f) a través dunha fórmula de cuadratura de tipo interpolatorio-polinómico In(f) de n+ 1 puntos, non se comete erro se a función f ten grao ≤n . O seguinte resultado proporciona a expresión do erro nas fórmulas de Newton-Cotes, no caso de que este exista. Denotarase πn(x) = n Y i=0 (x−xi)=(x−x0)(x−x1)... (x−xn), onde x0, x1, ... , xn son os n+ 1 nodos da fórmula de cuadratura. Teorema 1.9. Sexa unha función f∈ Cn+2([a, b]) e unha fórmula de Newton-Cotes con n par. O seu erro vén dado pola expresión En(f) = Mn (n+ 2)!hn+3f(n+2)(ξ), ξ ∈(a, b), (1.10) de modo que se a fórmula é pechada Mn=Zn 0 t πn(t)dt < 0, e se é aberta, Mn=Zn+1 −1 t πn(t)dt > 0. Sexa agora unha función f∈ Cn+1([a, b]) e unha fórmula de Newton-Cotes cun valor impar para n . O seu erro vén dado por En(f) = Kn (n+ 1)!hn+2f(n+1)(ξ), ξ ∈(a, b), (1.11) 1.2. Fórmulas de Newton-Cotes 9 de modo que se a fórmula é pechada Kn=Zn 0 πn(t)dt < 0, e se é aberta Kn=Zn+1 −1 πn(t)dt > 0. A demostración deste teorema pódese atopar en Isaacson e Keller (1994). A partir deste resultado pódese deducir o grao de exactitude das fórmulas de Newton-Cotes para os casos n par e impar, xa que non é o mesmo posto, que os seus erros de cuadratura non se calculan da mesma forma. Por ser as fórmulas de Newton-Cotes do tipo interpolatoriopolinómico, xa se viu no Teorema 1.5 que, para n+ 1 puntos, o seu grao de exactitude é polo menos n . Se n é impar, na expresión (1.11) vese que o erro é nulo cando f(n+1) = 0 , polo que a fórmula ten grao de exactitude n , como se acaba de deducir a partir do Teorema 1.5. En cambio, para o caso n par, razonando de forma análoga a partir de (1.10), dedúcese que a fórmula ten grao de exactitude n+ 1 , xa que, como o erro é nulo para f(n+2) = 0 , a fórmula de cuadratura é exacta para polinomios de grao ≤n+ 1 . 1.2.2. Obtención das fórmulas de Newton-Cotes A continuación móstranse os casos particulares de fórmulas de Newton-Cotes pechadas que se èmpregan máis habitualmente na práctica. Están aplicadas a un intervalo de integración [a, b] arbitrario, pechado e acotado. - Fórmula do trapecio (n=1) Neste caso tómanse dous nodos x0=a e x1=b , h=b−a , e os pesos de cuadratura son α0=α1=a+b 2 . A fórmula de cuadratura resultante é I1(f) = h 2[f(x0) + f(x1)], e, para f∈ C2([a, b]) , o de erro de cuadratura vén dado por E1(f) = −h3 12f00(ξ), ξ ∈(a, b). Esta fórmula ten grao de exactitude 1. - Fórmula de Simpson (n=2) Tómanse tres nodos x0=a , x1=a+b 2 e x2=b con h=b−a 2 , e os pesos α0=α2=b−a 6 e α1=2(b−a) 3 . A fórmula de cuadratura é I2(f) = h 3[f(x0)+4f(x1) + f(x2)], 10 1. Preliminares e, para f∈ C4([a, b]) , ten erro de cuadratura E2(f) = −h5 90f(4)(ξ), ξ ∈(a, b). O grao de exactitutde é 3. - Fórmula de Simpson 3/8 (n=3) Neste caso, tómanse catro nodos x0=a , x1=2a+b 3 , x2=a+2b 3 e x3=b , h= (b−a)/3 , e os pesos de cuadratura son α0=α3=b−a 8 , α1=α2=3(b−a) 8 , polo que a fórmula queda como I3(f) = 3h 8[f(x0)+3f(x1)+3f(x2) + f(x3)], que, para f∈ C4([a, b]) , ten erro de cuadratura E3(f) = −3h5 80 f(4)(ξ), ξ ∈(a, b). O grao de exactitude é 3. No seguinte exemplo, compáranse as fórmulas do trapecio e de Simpson que teñen dous e tres puntos, respectivamente, para a función f(x) = ex . Ademáis, tamén se consideran dous intervalos de integración de diferentes lonxitudes para mostrar as variacións do erro de cuadratura cando se varían estes dous factores. Exemplo 1.10. Comparación das fórmulas do trapecio e de Simpson para as integrais Z2 0 exdx e Z4 0 exdx : (a) O valor real da primeira integral é I(f) = Z2 0 exdx = 6.38906 e, para as fórmulas de cuadratura consideradas, obtense Trapecio: I1(f) = 1[f(x0) + f(x1)] = e0+e2= 8.38906 Simpson: I2(f) = 1 3[f(x0) + f(x1) + f(x2)] = 1 3(e0+ 4e1+e2)=6.42073 Pode verse que o erro se reduce considerablemente, pasando de −2 a −0.03167 , tomando únicamente un punto máis. (b) No segundo caso, o valor real da integral é I(f) = Z4 0 exdx = 53.59815 , e para as fórmulas de cuadratura consideradas, Trapecio: I1(f) = 2[f(x0) + f(x1)] = 2(e0+e4) = 111.19630 Simpson: I2(f) = 2 3[f(x0) + f(x1) + f(x2)] = 2 3(e0+ 4e2+e4) = 56.76958 Volve reducirse o erro de forma considerable, xa que pasa de −57.59815 a −3.17143 . Nesta ocasión, ao tomar un intervalo de integración maior, o erro de cuadratura obtido coa fórmula de Simpson segue sendo elevado, polo que sería favorable considerar unha fórmula de máis puntos.  1.2. Fórmulas de Newton-Cotes 11 No caso aberto, destácanse os seguintes casos particulares. - Fórmula do punto medio (n=0) Tense un único nodo, x0=b+a 2 , que se calcula como o punto medio do intervalo de integración [a, b] , h=b−a 2 e o peso é α0=b−a . Entón a fórmula resultante é I0(f)=2hf(x0), que para C2([a, b]) ten erro de cuadratura E0(f) = h3 3f00(ξ), ξ ∈(a, b). O seu grao de exactitude é 1. - n=1 Os nodos de cuadratura son x0=2a+b 3 , x1=a+2b 3 , h=b−a 3 , e os pesos α0=α1=b−a 2 , polo que, para C2([a, b]) , a fórmula é I(f) = 3h 2[f(x0) + f(x1)] + 3h3 4f00(ξ), ξ ∈(a, b) con grao de exactitude 2. - n=2 Os nodos son x0=3a+b 4 , x1=a+b 2 e x2=a+3b 4 , h=b−a 4 , e os pesos α0=α2=2(b−a) 3 e α1=b−a 3 . Polo tanto, para C4([a, b]) , I(f) = 4h 3[2f(x0)−f(x1)+2f(x2)] + 14h5 45 f(4)(ξ), ξ ∈(a, b). Ten grao de exactitude 3. No Exemplo 1.10 ilústrase un dos principais problemas que xorden cando se empregan as fórmulas de Newton-Cotes. Ao considerar intervalos de integración moi longos, tamén é necesario considerar moitos puntos no intervalo para unha correcta aproximación do valor real da integral. Pero xa se viu que cando o número de nodos n é elevado, os pesos αi , i= 0, ... , n poden tomar valores negativos, e que polo tanto un pequeno erro nas evaluacións f(xi) poida derivar en erros moi altos na aproximación da integral, causando así inestabilidades. Polo tanto, non por aumentar o número de nodos, se vai conseguir converxencia nestas fórmulas. É por este motivo polo que, alternativamente a estas, se empregan as fórmulas de Newton- Cotes compostas. A súa obtención é similar a das simples, pero non é necesario empregar fórmulas de graos moi altos para obter unha correcta aproximación da integral, o que evita as inestabilidades mencionadas anteriormente. 18 2. Fórmulas de Gauss sendo j, k = 0,1,2, ... . Dirase que esta sucesión é ortonormal se ademais para todo j, k = 0,1, ... se ten hpj, pkiw=Zb a w(x)pj(x)pk(x)dx = 1, para j =k. A proposición que se enuncia a continuación, vai ser moi importante para a determinación dos polinomios ortogonais. Proposición 2.5. Un polinomio pn(x) ortogonal respecto de w(x) , que debe satisfacer a Denición 2.3, existe é único salvo factor multiplicativo. Demostración. Verase que é posible atopar a0, ... , an−1 reais de modo que Zb a w(x)[a0+a1x+···+an−1xn−1+xn]xkdx = 0 para k= 0, ... , n −1 . Polo tanto téñense as seguintes n ecuacións (unha para cada valor de k ):                      Zb a w(x)[a0+a1x+···+an−1xn−1+xn]dx = 0 Zb a w(x)[a0x+a1x2+···+an−1xn+xn+1]dx = 0 . . . Zb a w(x)[a0xn−1+a1xn+···+an−1x2n−2+x2n−1]dx = 0. (2.2) Recórdese que, para cada k= 0,1, ... , 2n−1 , os momentos están denidos mediante µk=Zb a w(x)xkdx. Polo tanto, as ecuacións (2.2) pódense escribir como o seguinte sistema linear              a0µ0+a1µ1+... +an−1µn−1=−µn a0µ1+a1µ2+... +an−1µn=−µn+1 . . . a0µn−1+a1µn+... +an−1µ2n−2=−µ2n−1 (2.3) onde a0, ... , an−1 son as incógnitas. A continuación verase que o sistema (2.3) é non singular e que ten unha única solución. Para probar esto, suporase que o sistema (2.3) é singular, e que polo tanto  µ0µ1··· µn−1 µ1µ2··· µn . . .. . .. . . µn−1µn··· µ2n−2  = 0. (2.4) 2.1. Polinomios ortogonais 19 Isto equivale a atopar n constantes reais b0, b1, ... , bn−1 , non todas iguais a cero, tales que exista unha combinación linear das las de (2.4) igual a cero, e que polo tanto se teña o seguinte sistema linear              b0µ0+b1µ1+... +bn−1µn−1= 0 b0µ1+b1µ2+... +bn−1µn= 0 . . . b0µn−1+b1µn+... +bn−1µ2n−2= 0 ou equivalentemente,                      Zb a w(x)[b0+b1x+···+bn−1xn−1+xn]dx = 0 Zb a w(x)[b0x+b1x2+···+bn−1xn+xn+1]dx = 0 . . . Zb a w(x)[b0xn−1+b1xn+···+bn−1x2n−2+x2n−1]dx = 0. Multiplicando a primeira destas ecuacións por b0 , a segunda por b1 , e así sucesivamente, obtense que Zb a w(x)[b0+b1x+···+bn−1xn−1+xn]2dx = 0 (2.5) Como b0, b1, ... , bn−1 non son todas iguais a cero, o polinomio xn+bn−1xn−1+···+b1x+b0 é distinto de cero. Polo tanto, por como se viu na Denición 2.1 que debe ser a función peso w(x) , Zb a w(x)[xn+bn−1xn−1+···+b1x+b0]2dx > 0. Esto último contradí a expresión (2.5), polo que as constantes b0, b1, ... , bn−1 tal e como foron escollidas non poden existir. Polo tanto, o sistema de ecuacións (2.2) é non singular e ten unha única solución, polo que os coecientes a0, a1, ... , an−1 son únicos, o que implica que o polinomio xn+an−1xn−1+···+a1x+a0 tamén o sexa, salvo factor multiplicativo. Este resultado, implica que, sen perda de xeneralidade, se poidan considerar polinomios ortogonais pn(x) de grao n, n ≥1 , mónicos (é dicir, con coeciente principal igual a 1), da forma pn(x) = xn+an−1xn−1+···+a1x+a0, onde a0, ... , an−1 son números reais, e que ademais as sucesións de polinomios ortogonais tamén sexan únicas salvo factor multiplicativo. Polo tanto, a paartir de agora os polinomios ortogonais que se van considerar serán mónicos. 20 2. Fórmulas de Gauss A proba da Proposición 2.5 proporciona un método para o cálculo dos polinomios ortogonais para unha función peso e un intervalo de integración determinados. Ilústrase este proceso no seguinte exemplo. Exemplo 2.6. Cálculo do polinomio de Legendre de 3 º grao. Considéranse a función peso w(x)=1 e o intervalo de integración [−1,1] . Os momentos calcúlanse mediante µk=Zb a w(x)xkdx, k = 0,1, ... , 2n−1 e neste caso correspóndense con µ0= 2, µ2=2 3, µ4=2 5, µ1=µ3=µ5= 0. Polo tanto, o sistema de ecuacións correspondente (calculado de igual modo que (2.3)) é            2a0+ 0a1+2 3a2= 0 0a0+2 3a1+ 0a2=−2 5 2 3a0+ 0a1+2 5a2= 0 de onde se obteñen os valores a0= 0, a1=−3 5, a2= 0 para os coecientes de P(0) 3(x) . Chégase a que o polinomio buscado é P(0) 3(x) = x3−3 5x.  Cando os valores de n son altos, a resolución dun sistema como o do Exemplo 2.6 resulta un proceso moito máis complexo. O seguinte resultado, que establece unha relación de recurrencia para o cálculo dos polinomios ortogonais, fai que estes se poidan calcular dunha forma moito máis sinxela. Ademais, posteriormente tamén se verá que esta relación será de gran importancia para a obtención dos pesos e nodos de cuadratura nas fórmulas gaussianas. Teorema 2.7. Sexa unha función peso w(x) denida en [a, b] , a sucesión de polinomios ortogonais {pn}∞ n=0 satisfai a seguinte relación de recurrencia de tres termos: pk+1(x)=(x−ak)pk(x)−bkpk−1(x)k≥0, p−1(x) = 0, p0(x)=1, (2.6) 2.1. Polinomios ortogonais 21 de modo que ak=hxpk, pkiw hpk, pkiw k≥0, (2.7) bk=hpk, pkiw hpk−1, pk−1iw k≥1, (2.8) denotando b0=Zb a w(x)dx . Demostración. Os polinomios construiranse mediante a ortogonalización de Gram-Schmidt,e seguirase un proceso inductivo partindo de que p0(x) = 1 . Supóñase que, por hípotese de indución, todos os polinomios ortogonais como os anteriores están contruídos ata certo i≤k . Pola Proposición 2.5, estes polinomios son únicos para cada i= 0, ... , k . Queda probar que existe un único polinomio de grao k+ 1 tal que hpk+1, piiw= 0 para todo i≤k, (2.9) que satisfaga a relación de recurrencia (2.6). Calquera polinomio pk+1 ∈Pk+1 pódese escribir, dadas c0, ... , ck−1 constantes reais, como pk+1(x)=(x−ak)pk(x) + ck−1pk−1(x) + ck−2pk−2(x) + ... +c0p0(x), xa que o seu coeciente principal e os de todos os polinomios pi, i ≤k son iguais a 1, posto que sa se viu que son mómicos. Como ademais hpi, pjiw= 0 para j≤k, i 6=j , partindo de (2.9), obtense que hpk+1, pkiw=hpk, xpkiw−akhpk, pkiw= 0 (2.10) hpk+1, pi−1iw=hxpi−1, pkiw+ci−1hpi−1, pi−1iw= 0, i ≤k. (2.11) Téñense os produtos escalares hpk, pkiw= 0 e hpi−1, pi−1iw= 0 para i= 1, ... , k , xa que p2 k e p2 i−1 son polinomios non negativos. Polo tanto, as dúas ecuacións anteriores teñen unha única solución. Despexando na primeira delas obtense a expresión (2.7) para os coecientes ak . Agora, pola hipótese de indución, da relación de recurrencia (2.6), pi(x)=(x−ai−1)pi−1(x)−bi−1pi−2(x), i ≤k. Polo tanto, resolvendo esto para xpi−1(x) , obtense que hxpi−1, pkiw=hpi, pkiw para i≤k , e despexando en (2.11), ci−1=−hpi, pkiw hpi−1, pk−1iw =(−bkse i =k 0se j < k, e obtense a expresión (2.8) para o cálculo dos coecentes bk . 22 2. Fórmulas de Gauss Os coecientes ak e bk que se acaban de denir no Teorema 2.7 son constantes e reais, con bk>0 para todo k= 0, ... , n , e permiten representar a sucesión de polinomios de forma única e compacta. A seguinte proposición, é unha propiedade importante sobre os ceros dos polinomios ortogonais, xa que estes ceros van ser os nodos na fórmula de cuadratura gaussiana correspondente. Proposición 2.8. Sexa pn(x) o polinomio de grao n ortogonal con todos os polinomios de grao ≤n−1 no intervalo [a, b] para w(x) . Entón, todos os ceros de pn(x) son reais, simples, e están no intervalo aberto (a, b) . Demostración. Supóñase que pn(x) ten un total de s≤n ceros diferentes, de modo que se pode escribir este polinomio como pn(x)=(x−x0)m0(x−x1)m1... (x−xs−1)ms−1, onde mi é a multiplicidade de xi para i= 0, ... , s −1 , e m0+... +ms−1=n . Ademais, tamén se suporá que estes ceros están nomeados de modo que x0, ... , xk−1 , k≤s son os ceros con multiplicidade impar, é dicir, aqueles nos que pn(x) cambia de signo. Se k=n , todos os ceros de pn(x) son reais, distintos (é dicir, simples), e están en (a, b) . Polo tanto queda probada a proposición. Se pola contra, se supón que k < n , verase a continuación que non pode ser posible. Neste caso, tómase o polinomio qk(x)=(x−x0)(x−x1)... (x−xk−1) (se k= 0 , considerarase q0(x)=1 ), que cambia de signo en (a, b) nos mesmos puntos que pn(x) , polo que o produto qk(x)pn(x) non cambia de signo en (a, b) . Entón, Zb a w(x)qk(x)pn(x)dx 6= 0. (2.12) Pero, como k < n , e pn(x) é ortogonal con todos os polinomios de grao ≤n−1 , tense que, por ser qk(x) un polinomio de grao k≤n−1 , Zb a w(x)qk(x)pn(x)dx = 0. Isto contradí a expresión (2.12), e polo tanto, k < n non pode ser posible, e necesariamente os n ceros de pn(x) son simples. A continuación móstranse as principais familias clásicas polinomios ortogonais, que están asociadas ás funcións peso e intervalos do Cadro 2.1, e para os cales os coecientes ak e bk son coñecidos de forma explícita, e polo tanto, tamén a sucesión de polinomios ortogonais asociada. 2.1. Polinomios ortogonais 23 Polinomios ortogonais de Jacobi Recórdese que a súa función peso asociada é w(x) = (1−x)α(1+x)β, α, β > −1 , e o intervalo de integración é [−1,1] . Os polinomios de Jacobi de grao n adoitan denotarse por P(α,β) n(x) . Pola Proposición 2.8, sábese que os seus n ceros son reais e simples, e que están situados sobre o intervalo aberto (−1,1) . Variando os valores de α e de β , obtéñense outras dúas familias de polinomios, que son as máis usadas nas cuadraturas gaussianas: Os polinomios de Legendre danse cando α=β= 0 , polo que a función peso é w(x)=1 e o intervalo de integración, [−1,1] . Denotaranse por P(0) n(x) . A relación de recurrencia é neste caso P(0) k+1(x) = xP (0) k(x)−k2 (2k−1)(2k+ 1)P(0) k−1(x), k ≥0 polo que os correspondentes coecientes son ak= 0, k ≥0 e bk=k2 (2k−1)(2k+ 1),1k≥1. (2.13) Os polinomios de Chebyshev de primeira especie obtéñense tomando α=β=−1/2 . Adoitan denotarse por Tn(x) . Os coecientes da relación de recurrencia son ak= 0 para k≥0 , e bk=(1/2, se k= 1 1/4, se k≥2. Polinomios ortogonais de Laguerre Neste caso, a función peso é w(x) = xαe−x, α > −1 e o intervalo [0,∞) . Estes polinomios denótanse por L(α) n(x) . Pola Proposición 2.8, as súas n raíces son reais e simples están situadas no intervalo aberto (0,∞) . Neste caso a relación de recurrencia que se establece é L(α) k+1(x) = (2k+α+ 1 −x)L(α) k(x)−k(k+α)L(α) k−1(x), k ≥0, e polo tanto, os coecientes son ak= 2k+α+ 1, k ≥0 e bk=k(k+α), k ≥1 . Polinomios ortogonais de Hermite Considéranse a función peso w(x) = e−x2 e o intervalo (−∞,∞) . Os polinomios de Hermite de grao n denótanse por Hn(x) . Pola Proposición 2.8, os seus n ceros son reais e simples e sitúanse en todo R , xa que é o intervalo de integración. Polo tanto, a relación de recurrencia queda dada por Hk+1(x) = xHk(x)−k 2Hk−1(x), k ≥0, onde os coecientes son ak= 0, k ≥0 e bk=k 2, k ≥1 . 24 2. Fórmulas de Gauss 2.2. Fórmulas de cuadratura de Gauss 2.2.1. Máximo grao de exactitude A principal vantaxe á hora de usar as fórmulas de cuadratura gaussiana é que se alcanza un grao de exactitude de 2n+ 1 , que, como se verá, é o máximo posible cando se está a falar de fórmulas integración numérica de n+ 1 puntos. Agora mostrarase que as fórmulas de cuadratura construídas cos nodos e os pesos calculados como se viu na sección anterior alcanzan este grao máximo de exactitude, e que non é posible superalo. Para probar esto danse os seguintes resultados. Antes, recórdase a denición do seguinte polinomio asociado aos n+ 1 nodos de cuadratura x0, ... , xn , xa que se empregará de novo en algunhas demostracións: πn= n Y i=0 (x−xi). Teorema 2.9. Sexa k un número natural tal que 1≤k≤n+ 1 . Unha fórmula de cuadratura da forma I(f) = Zb a w(x)f(x)dx = n X i=0 αif(xi) + En(f) ten grao de exactitude polo menos n+k se e só se se cumpren as seguintes condicións: 1. A fórmula de cuadratura é de tipo interpolatorio-polinómico. 2. Para todo polinomio q(x)∈Pk−1 se cumpre hπn, qiw=Zb a w(x)πn(x)q(x)dx = 0, (2.14) é dicir que πn(x) e q(x) son polinomios ortogonais respecto de w(x) . Demostración. (⇒) A primeira condición próbase no Teorema 1.5. Para probar a segunda condición, considerarase un polinomio q(x)∈Pk . Polo tanto, o máximo grao que poderá alcanzar o polinomio πn(x)·q(x) será n+k . Por hipótese, a fórmula de cuadratura ten grao de exactitude n+k , logo tense Zb a w(x)πn(x)q(x)dx = n X i=0 αiπn(xi)q(xi)=0, que vale 0 porque πn(xi) = 0 para i= 0, ... , n . Polo tanto, pola denición de ortogonalidade, pode dicirse que πn(x) é ortogonal respecto de w(x) con todos os polinomios de grao ≤n . 2.2. Fórmulas de cuadratura de Gauss 25 (⇐) Sexa p(x) un polinomio de grao ≤n+k para certo valor 1≤k≤n+ 1 . Dividindo este polinomio p(x) entre πn(x) , p(x) = q(x)πn(x) + r(x), q ∈Pk−1, r ∈Pn. Tendo en conta que πn(xi) = 0 para todo i= 0, ... , n , tense que p(xi) = q(xi)πn(xi) + r(xi) = r(xi), i = 0, ... , n. Usando a condición 2, é dicir, que πn(x) é ortogonal con todos os polinomios de grao menor ou igual que n , Zb a w(x)πn(x)q(x)dx = 0. (2.15) Como a fórmula de cuadratura é de tipo interpolatorio-polinómico e ten n+ 1 puntos, xa se viu no primeiro capítulo, no Teorema 1.5, que é exacta para polinomios de grao n , logo Zb a w(x)r(x)dx = n X i=0 αir(xi). (2.16) Polo tanto, a partir destas últimas expresións (2.15) e (2.16), tense que Zb a w(x)p(x)dx =Zb a w(x)πn(x)q(x)dx +Zb a w(x)r(x)dx = =Zb a w(x)r(x)dx = n X i=0 αir(xi) = n X i=0 αip(xi), e pódese ver que a fórmula de cuadratura é exacta para calquera polinomio de grao n+k , onde 1≤k≤n+ 1 . Corolario 2.10. Nas condicións do Teorema 2.9, o máximo grao de exactitude da fórmula de cuadratura é 2n+ 1 , e conséguese se e so se os n+ 1 nodos son os ceros do polinomio ortogonal de grao n+ 1 respecto da función peso w(x) en [a, b] e os pesos de cuadratura verican αi=Zb a w(x)li(x)dx, i = 0, ... , n, onde li(x) son os polinomios característicos de Lagrange, que xa se deniron en (1.3) . Demostración. Tomando k=n+ 1 , tense que o grao de exactitude da fórmula é 2n+ 1 . Verase que existe polo menos un polinomio de grao 2n+ 2 para o cal a fórmula non é exacta. Considerarase o polinomio πn(x)∈Pn+1 , polo tanto π2 n(x)∈P2n+2 . Claramente, por como xa se viu que debe ser a función peso w(x) , Zb a w(x)π2 n(x)dx > 0, (2.17) 26 2. Fórmulas de Gauss xa que π2 n(x) é un polinomio positivo en [a, b] . Pero por outra parte, estase supoñendo que a fórmula de cuadratura é exacta para π2 n(x) , logo Zb a w(x)π2 n(x)dx = n X i=0 αiπ2 n(xi)=0, (2.18) o cal é unha contradición con (2.17). Polo tanto a fórmula de cuadratura non pode ter grao de exactitude 2n+ 1 . Nótese que na expresión (2.18), a segunda igualdade dáse porque x0, ... , xn son os n+1 ceros do polinomio πn(x) . Tamén en relación coa converxencia das fórmulas de Gauss, enúnciase o seguinte resultado sobre os pesos de cuadratura, que será importante ter en conta cando se fale da converxencia desta familia de fórmulas. Proposición 2.11. Nas condicións do Teorema 2.9, todos os pesos de cuadratura nas fórmulas de Gauss son positivos. Ademais, cúmprese α0+α1+... +αn=Zb a w(x)dx. A proba da primeira parte pódese atopar en Stroud (1974), e da segunda en Martínez- Finkelshtein e Moreno-Balcázar (1999). Recórdese que nas fórmulas de Newton-Cotes, tanto nas pechadas como nas abertas, a partir de certos valores de n existía a presenza de pesos de cuadratura negativos, e que este feito ocasionaba irregularidades en canto converxencia do método. Coas fórmulas de Gauss, isto non sucede, posto que a Proposición 2.11 asegura a positividade de todos os pesos de cuadratura. 2.2.2. Unicidade das fórmulas de Gauss Agora probarase a unicidade das fórmulas de Gauss, é dicir, que para un intervalo [a, b] e unha función peso w(x) determinadas, existe unha única fórmula de cuadratura de n+ 1 puntos con grao de exactitude 2n+ 1 . Teorema 2.12. Se unha fórmula de cuadratura é exacta para todos os polinomios de grao ≤ 2n+ 1 , os nodos deben ser os ceros do polinomio pn+1(x) , ortogonal con todos os polinomios de grao ≤n para o intervalo de integración [a, b] e a función peso w(x) . Demostración. Sexa unha fórmula In(f) = n X i=0 αif(xi), (2.19) 2.3. Cálculo das fórmulas de cuadratura gaussianas 27 exacta para polinomios de grao ≤2n+ 1 . Considerando o polinomio πn(x) = n Y i=0 (x−xi), se qn(x) é un polinomio arbitrario de grao n , a fórmula (2.19) será exacta para πn(x)qn(x) . Polo tanto, como πn(x) = 0 para i= 0, ... , n , tense que Zb a w(x)πn(x)qn(x)dx = n X i=0 αiπn(xi)qn(xi)=0. Equivalentemente a esta última expresión, pode dicirse que πn(x) é ortogonal con todos os polinomios de grao n . Pola Proposición 2.5 da unicidade dos polinomios ortogonais, dado un intervalo de integración e unha función peso determinados, tense que πn(x) = pn+1(x) , onde pn+1(x) é o polinomio ortogonal de grao n+ 1 asociado ao intervalo [a, b] e á función peso w(x) . 2.2.3. Estimación do erro nas fórmulas de Gauss O teorema que se enuncia a continuación, cuxa proba se pode atopar en Stoer e Bulirsch (2010), dá a estimación do erro para as cuadraturas gaussianas. Teorema 2.13. Sexa unha función f∈ C2n+2([a, b]) , o erro de cuadratura da fórmula está dado por En(f) = Zb a w(x)f(x)dx − n X i=0 αif(xi) = f(2n+2)(ξ) (2n+ 2)! hpn+1, pn+1iw, ξ ∈(a, b). (2.20) Neste tipo de fórmulas, a diferenza das de Newton-Cotes, canto maior é o número de nodos, menor é o erro cometido. É dicir, En(f)−→ 0 cando n−→ ∞ . Isto implica que os resultados obtidos coas fórmulas de Gauss sexan máis exactos que os de Newton-Cotes, xa que o grao de precisión se duplica. Este feito unido a que, como xa se viu na Proposición 2.11, os pesos de cuadratura das fórmulas gaussianas son sempre positivos, fai que as aproximacións obtidas sexan considerablemente mellores. O problema que se atopa coas fórmulas de Gauss tratarase a continuación, e é a dicultade no cálculo dos pesos e dos nodos de cuadratura. 2.3. Cálculo das fórmulas de cuadratura gaussianas Ata este momento, viuse que para obter as fórmulas de cuadratura gaussianas é necesario, primeiro calcular o polinomio ortogonal asociado ao intervalo e á función peso que correspondan, e a partir deste calcular os nodos e os pesos da fórmula. Os nodos recórdese que se corresponden 34 3. Integración de Romberg 3.1. Extrapolación de Richardson A extrapolación de Richardson é un método de derivación numérica a través do cal se pretende mellorar o cálculo da derivada dunha función determinada diminuíndo o erro que se obtén a través do parámetro h . Para isto, combínanse diversas aproximacións de baixa orde de dita función para obter unha mellor aproximación con maior orde. Considérese, por exemplo, unha función N(h) , denida por N(h) = N0+khm+O(hm+1), (3.1) onde N(h) é unha aproximación de N0 de orde m e h tende a 0 . O propósito da extrapolación de Richardson é obter unha nova aproximación da cantidade N0 de maior orde. Neste caso será de orde m+ 1 . Substituíndo agora o parámetro h por δh para certo 0< δ < 1 , reescríbese a expresión (3.1) como N(δh) = N0+k(δh)m+O(hm+1). (3.2) Multiplicando (3.1) por δm e restándolle (3.2), obtense δmN(h)−N(δh)=(δm−1)N0+ (δm−1)O(hm+1) δmN(h)−N(δh) δm−1=N0+O(hm+1). Obsérvese que empregando este método, obtéñense aproximacións da orde que se desexe, xa que continuando o proceso que se acaba de explicar, pódense ir desarrollando métodos de orde O(hm+2),O(hm+3), ... , ata obter unha aproximación para N0 tan precisa como se desexe. 3.2. Integración de Romberg Este método de cuadratura numérica está baseado na extrapolación de Richardson. Recórdese que este algoritmo se aplica para fórmulas de derivación numérica. Nesta ocasión a idea é calcular diversas aproximacións preliminares dunha integral I(f) = Zb a f(x)dx mediante a fórmula do trapecio composta para distintos valores de n . Despois, mediante o proceso de extrapolación de Richardson, mellóranse estas primeiras aproximacións para obter un mellor resultado, cuxa orde de converxencia sexa maior. Xa se estudou no capítulo anterior que a fórmula do trapecio composta ten un erro de orde O(h2) , e que, para unha función f∈ C∞([a, b]) , xando un número m de subintervalos de [a, b] , 3.2. Integración de Romberg 35 m≥1 , e tomando h=b−a m e os nodos xj=a+jh para j= 0, ... , m , ten a forma Zb a f(x)dx =h 2 f(a)+2 m−1 X j=1 f(xj) + f(b) −b−a 12 h2f00(ξ), ξ ∈(a, b). Como xa se mencionou, para levar a cabo este proceso, primeiro é necesario calcular r aproximacións preliminares da integral I(f) mediante a fórmula do trapecio composta, nas que o número de subintervalos vén dado pola cantidade m= 2j , para j= 0, ... , r . Polo tanto, hj=b−a 2j e os nodos son xj=a+jh para cada j= 0, ... , r . Deste modo, as r+ 1 fórmulas do trapecio necesarias para iniciar a integración de Romberg están dadas por I0=h0 2[f(a) + f(b)] −b−a 12 h2 0f00(ξ0), ξ0∈(a, b) I1=h1 2[f(a)+2f(a+b 2) + f(b)] −b−a 12 h2 1f00(ξ1), ξ1∈(a, b) . . . Ir=hr 2 f(a)+2 2r−1 X j=1 f(xj) + f(b) −b−a 12 h2 rf00(ξr), ξr∈(a, b). Tendo en conta que Ij=R(j, 0) + Ej para j= 0, ... , r , onde se denota Ij como a fórmula do trapecio composta de 2j subintervalos de [a, b] , R(j, 0) como a aproximación obtida no proceso e Ej como o erro de cuadratura cometido nela, calcúlanse os valores de R(j, 1) como segue. Este será o primeiro paso do proceso. Para un valor j= 1, ... , r arbitrario considéranse Ij−1=R(j−1,0) + Ej−1= =hj−1 2 f(a)+2 2j−1−1 X i=1 f(xi) + f(b) −b−a 12 h2 j−1f00(ξj−1), ξj−1∈(a, b) Ij=R(j, 0) + Ej= =hj 2 f(a)+2 2j−1 X i=1 f(xi) + f(b) −b−a 12 h2 jf00(ξj), ξj∈(a, b). Tomando a derivada f00 como unha función constante, tense que f00(ξj−1) = f00(ξj) , e polo tanto Ej−1 Ej=−b−a 12 h2 j−1f00(ξj−1) −b−a 12 h2 jf00(ξj) =hj−1 hj2 ⇒Ej−1=Ejhj−1 hj2 . 36 3. Integración de Romberg Tendo en conta que I(f) = Ij=Ij−1⇒R(j, 0) + Ej=R(j−1,0) + Ej−1 , R(j, 0) + Ej=R(j−1,0) + Ejhj−1 hj2 ⇒Ej=R(j−1,0) −R(j, 0) 1−hj−1 hj2. (3.3) Ademais, é trivial que hj−1= 2hj . Logo substituíndo esto na expresión (3.3), obtense Ej=R(j−1,0) −R(j, 0) 1−22=R(j, 0) −R(j−1,0) 3 R(j, 0) + Ej=R(j, 0) + R(j−1,0) −R(j, 0) 1−22=Ij+Ij−Ij−1 3=4 3R(j, 0) + 1 3R(j−1,0). Así, para j= 1, ... , r , tense R(j, 1) = 4 3R(j, 0) −1 3R(j−1,0), que son r novas aproximacións da integral I(f) de orde O(h4) . O seguinte paso do proceso consiste na obtención das r−1 aproximacións R(j, 2) para j= 2, ... , r . Débese ter en conta que que se calcularon previamente teñen orde de converxencia O(h4) , polo tanto o seu erro terá a forma Ch4 , onde C é unha constante real. Seguindo un procedemento análogo ao anterior, obtéñense no seguinte paso da extrapolación as r−1 aproximacións de orde O(h6) , que se calculan mediante R(j, 2) = 16 15R(j, 1) −1 15R(j−1,1), j = 2, ... , r. En xeral as aproximacións de I(f) están dadas por R(j, k) = 1 4k−1(4kR(j, k −1) −R(j−1, k −1)), j = 1, ... , r k=j, ... , r e teñen orde O(h2(k+1)) , polo que no último paso deste proceso, obterase unha aproximación R(r, r) con orde de converxencia O(h2(r+1)) . Na seguinte táboa ilústranse os pasos seguidos neste proceso, ademais da orde de converxencia que se obtén en cada un dos pasos da extrapolación. Trapecio Extrapolación Extrapolación Extrapolación jO(h2)O(h4)O(h6)··· O(h2(r+1)) 0 R(0,0) 1 R(1,0) R(1,1) 2 R(2,0) R(2,1) R(2,2) . . .. . .. . .. . .... r R(r, 0) R(r, 1) R(r, 2) ··· R(r, r) 3.2. Integración de Romberg 37 Móstrase a continuación un exemplo que fai uso deste método que se acaba de explicar. Exemplo 3.2. Cálculo da integral I(f) = Zπ 0 sin (x)dx mediante o proceso de integración de Romberg, para r= 3 . Primeiro débense calcular as aproximacións preliminares mediante as fórmulas do trapecio compostas de 1, 2, 4 e 8 intervalos para a función f(x) = sin (x) no intervalo [0, π] : R(0,0) = π 2[ sin 0 + sin π] = 0 R(1,0) = π 4h sin 0+2 sin π 2+ sin πi= 1.57079633 R(2,0) = π 8 sin 0+2 sin π 4+ sin π 2+ sin 3π 4+ sin π= 1.89611890 R(3,0) = π 16  sin 0+2 sin π 8+···+ sin 7π 8+ sin π= 1.97423160 todas con orde de converxencia O(h2) . A continuación, empregando a fórmula R(j, 1) = 4 3R(j, 0) −1 3R(j−1,0), j = 1,2,3 calcúlanse as seguintes aproximacións, que terán orde de converxencia O(h4) : R(1,1) = 4 3R(1,0) −1 3R(0,0) = 2.09439511 R(2,1) = 4 3R(2,0) −1 3R(1,0) = 2.00455976 R(3,1) = 4 3R(3,0) −1 3R(2,0) = 2.00026917 Análogamente, para a obtención das aproximacións de orde O(h6) emprégase a fórmula correspondente R(j, 2) = 16 15R(j, 1) −1 15R(j−1,1), j = 2,3 e obtense R(2,2) = 16 15R(2,1) −1 15R(1,1) = 1.99857073 R(3,2) = 16 15R(3,1) −1 15R(2,1) = 1.99998313. Finalmente, a última aproximación, que terá orde de converxencia O(h8) , está dada por R(3,3) = 64 63R(3,2) −1 63R(2,2) = 2.00000555. Todo este proceso pódese resumir na seguinte táboa 38 3. Integración de Romberg jO(h2)O(h4)O(h6)O(h8) 0 0 1 1.57079633 2.09439511 2 1.89611890 2.00455976 1.99857073 3 1.97423160 2.00026917 1.99998313 2.00000555 tomando o valor R(3,3) = 2.00000555 como o resultado obtido na aproximación. Neste caso coñécese o valor real da integral, I(f) = Zπ 0 sin (x)dx = 2 , polo que se pode comprobar fácilmente que o erro cometido nesta aproximación é I(f)−R(3,3) = −0.00000555 . En caso de que se desexe aumentar a orde de converxencia do método, aumentarase o valor de r . Se por exemplo se toma r= 4 , calcúlanse os valores de R(4, k) para k= 1,2,3,4 , e táboa anterior queda como jO(h2)O(h4)O(h6)O(h8)O(h10) 0 0 1 1.57079633 2.09439511 2 1.89611890 2.00455976 1.99857073 3 1.97423160 2.00026917 1.99998313 2.00000555 4 1.99357034 2.00001659. 1.99999975 2.00000001 1.99999999 de modo que a aproximación neste caso é R(4,4) = 1.99999999 , con orde de converxencia O(h10) , e cun erro de I(f)−R(4,4) = 0.00000001 .  Bibliografía [1] Burden, R. L., Faires, D. J. e Burden, A. M. (2015). Numerical Analysis , 10th ed., Cengage Learning. [2] Conte, S. D. e de Boor, C. (1972). Analisis numerico elemental: Un enfoque algoritmico , 2 ª ed., McGraw-Hill. [3] Davis, P. J. e Rabinowitz, P. (1975). Methods of numerical integration , Computer Science ans Applied Mathematics. Academic Press, New York. [4] Flores, W. M. (2010). Introducción a los métodos numéricos . Instituto Tecnológico de Costa Rica. [5] Gautschi, W. (1994). Algorithm 726: ORTHPOLa package of routines for generating orthogonal polynomials and Gauss-type quadrature rules . ACM Transactions on Mathematical Software, 20 , 2162. [6] Gautschi W. (1982). On generating orthogonal polynomials . SIAM Journal on Scientic Computing, 3 , 289317. [7] Golub, G. H. e Welsch, J. H. (1969). Calculation of Gauss quadrature rules . Mathematics of Computation, 23 , 221230. [8] Isaacson, E. e Keller, H. B. (1994). Analysis of Numerical Methods . Dover Publications. [9] Martínez, J. J. (2001). Polinomios ortogonales, cuadratura gaussiana y problemas de valores propios . Margarita mathematica, 595606. [10] Martínez-Finkelshtein, A. e Moreno-Balcázar, J. J. (1999). Métodos numéricos: Aproximación en R . Almería, España: Serivicio de Publicaciones de la Universidad de Almería. [11] Mañas-Mañas J. F. e Pinta M. A. (2018). Métodos Numéricos para el Análisis Matemático con Matlab . UTMACH. 39 40 BIBLIOGRAFÍA [12] Quarteroni, A., Sacco, R. e Saleri, F. (2000). Numerical Mathematics , Texts in Applied Mathematics, 37. Springer-Verlag, New York. [13] Sanz-Senra, J.M. (2010). Diez lecciones de cálculo numérico , 2 ª ed.. Universidad de Valladolid. [14] Stoer, J. e Bulirsch, R. (2010). Introduction to Numerical Analysis , 3rd ed., Texts in Applied Mathematics, 12, Springer. [15] Stroud, A. H. (1974). Numerical Quadrature and Solution of Ordinary Dierential Equations , Applied Mathematical Sciences, 10. Springer-Verlag, New York. [16] Viaño, J. M., Burguera, M. (2000). Interpolación , Lecciones de métodos numéricos, 3, Tórculo Edicións. [17] Integración numérica . (s.f.). Facultad de Ingeniería de la Universidad Nacional del Nordeste. http://www.ing.unne.edu.ar/assets/pdf/academica/departamentos/computacion/comp/IN.pdf