Repositorio Institucional de Documentos

Trabajo Fin de Máster Factorización de matrices Vandermonde y fórmula de Newton Autora: Yasmina Khiar Viana Directores: Jesús Carnicer Álvarez Juan Manuel Peña Ferrández Máster en Modelización e Investigación Matemática, Estadística y Computación Zaragoza, Julio 2014 Facultad de Ciencias 2 Índice 1. Introducción 4 2. Problemas de interpolación y matrices de colocación 5 3. Interpolación polinómica y matriz de Vandermonde 6 4. Fórmula de Newton y factorización LU 7 5. Condicionamiento 10 6. Cálculo de las inversas de las matrices triangulares 10 7. Alta precisión relativa 12 8. La ordenación de Leja 13 9. Ejemplos 15 9.1. K= [0,1] ...................................... 15 9.2. K= [−1,1] ..................................... 17 10. Experimentos numéricos: factores triangulares 21 10.1. Intervalo [0,1] ................................... 21 10.2. Intervalo [−1,1] .................................. 24 11. Experimentos numéricos: condicionamiento conjunto 27 11.1. Intervalo [0,1] ................................... 28 11.2. Intervalo [−1,1] .................................. 29 12. Conclusiones 31 Bibliografía 33 3 1. Introducción Uno de los métodos más utilizados para resolver un sistema lineal de ecuaciones Ax=b consiste en descomponer la matriz Acomo producto de una matriz triangular inferior L y una matriz triangular superior Uy resolver los correspondientes sistemas triangulares. Esta forma de descomponer una matriz se conoce como factorización LU. En este trabajo vamos a considerar la matriz de Vandermonde Ven nodos x0, . . . , xn. Aunque disponemos de expresiones explícitas para la inversa de una matriz de Vandermonde, es habitual resolver el sistema mediante una factorización LU. Se explorarán las conexiones de la fórmula de interpolación de Newton con dichas factorizaciones triangulares. Como la fórmula de Newton depende de la ordenación de los nodos, se analizará el condicionamiento de factorizaciones triangulares asociadas a diferentes ordenaciones de los nodos. En cualquier método numérico, además del coste computacional, es importante tener en cuenta los errores de redondeo que se cometen. Una medida apropiada del error cometido es el condicionamiento de una matriz A, que se denota por κ∞(A)y se define del siguiente modo: κ∞(A) := kAk∞kA−1k∞. Es bien conocido que la matriz de Vandermonde está mal condicionada. Además, se tiene que κ∞(V)≤κ∞(L)κ∞(U). El objetivo de este trabajo es analizar posibles factorizaciones LU de Vy ordenaciones de los nodos de forma que κ∞(L)κ∞(U)sea lo menor posible. Las factorizaciones que vamos a considerar son dos: la primera está asociada a la fórmula de Newton para la resolución del correspondiente problema de interpolación polinómica y la segunda asociada a la eliminación gaussiana. Las ordenaciones que vamos a tomar son la natural, la de Leja y, en el caso de tomar los nodos en el intervalo [−1,1], proponemos una nueva ordenación llamada central. En la Sección 2 de este trabajo se introduce el problema de interpolación de Lagrange y se define la matriz de colocación de un problema de interpolación respecto de una base dada. A continuación en la Sección 3, se presenta el problema de interpolación de Lagrange polinómico y se da la solución de éste a través de la fórmula de Lagrange y de la matriz de Vandermonde. La Sección 4 comienza con la fórmula de Newton, que expresa el interpolante mediante diferencias divididas y se presenta la factorización LU asociada. Se dan fórmulas de recurrencia para calcular los elementos de las matrices triangulares LyU. Estas matrices triangulares forman una factorización LU de la matriz de Vandermonde Vdonde Utiene unos en la diagonal principal. En la siguiente Sección 5, se define el condicionamiento tradicional κ∞(A)y el condicionamiento de Skeel de una matriz. La Sección 6 está dedicada a calcular las inversas de las matrices LyUanteriores y se dan fórmulas explícitas que permiten calcular los elementos de estas matrices por recurrencia. En la siguiente Sección 7, se define el concepto de alta precisión relativa y se demuestra que los cálculos de L,L−1,UyU−1se realizan con alta precisión relativa. En la Sección 8 se introduce una nueva forma de ordenar los nodos: 4 la ordenación de Leja, descrita por Reichel en [10]. La Sección 9 está dedicada a ejemplos ilustrativos en los que se consideran distintos intervalos para tomar los nodos. Se consideran nodos equidistantes y se toman distintas ordenaciones para ellos. Las secciones 10 y 11 están dedicadas a los experimentos numéricos realizados. La primera de estas secciones recoge resultados sobre los factores triangulares, mientras que la segunda compara los resultados de los condicionamientos conjuntos. La última sección recoge las conclusiones extraídas de este trabajo. 2. Problemas de interpolación y matrices de colocación En este trabajo queremos estudiar la factorización de matrices Vandermonde. Estas matrices aparecen al plantear un problema de interpolación de Lagrange respecto de la base de monomios. Problema de interpolación de Lagrange. Sea Uun espacio de funciones con dim U= n+ 1. Dados nodos x0, . . . , xndistintos y f0, . . . , fn, encontrar u∈U, con u(xi) = fi, i= 0, . . . , n. El interpolante puede obtenerse en distintos espacios y dentro del mismo espacio puede expresarse en términos de diferentes bases. Si (u0, . . . , un)es una base ordenada de U, podemos expresar la solución uen términos de dicha base u=Pn i=0 aiuiy el problema de interpolación queda reducido a la resolución del sistema Mu0, . . . , un x0, . . . , xna=f, donde a= (a0, . . . , an)T,f= (f0, . . . , fn)Ty Mu0, . . . , un x0, . . . , xn= (uj(xi))i,j=0,...,n, denota la matriz de colocación de la base en los nodos. Aquí nos apartamos del convenio usual y consideramos que el primer índice de filas y columnas es 0en lugar de 1, por lo que la dimensión correspondiente a nes n+ 1, para acomodar nuestros resultados al problema de interpolación de grado ny a la notación de matrices de Vandermonde. Más generalmente podemos plantear problemas de interpolación lineales respecto a una sucesión de funcionales λ0, . . . , λn. Problema de interpolación lineal general. Sea Uun espacio de funciones de dimensión n+ 1. Dados funcionales λ0, . . . , λnyfuna función, encontrar u∈U, con λiu=λif, i= 0, . . . , n. 5 Introducimos la notación de matriz de colocación de un problema de interpolación respecto a una sucesión de funcionales. Definición. La matriz de colocación de una base (u0, . . . , un)respecto a una sucesión de funcionales (λ0, . . . , λn)es Mu0, . . . , un λ0, . . . , λn= (λi(uj))i,j=0,...,n. Una vez elegida una base, (u0, . . . , un), el problema de interpolación lineal general se reduce a resolver el sistema Mu0, . . . , un λ0, . . . , λna=Mf λ0, . . . , λn, expresando la solución en la forma u=Pn i=0 aiui. Un problema de interpolación lineal general no tiene necesariamente solución única. La condición necesaria y suficiente para que el problema admita una única solución es que el determinante de la matriz de colocación respecto de una base sea distinto de cero. 3. Interpolación polinómica y matriz de Vandermonde Comencemos esta sección planteando el problema de interpolación de Lagrange polinómico. Problema de interpolación de Lagrange polinómico. Dados nodos x0, . . . , xndistintos y f0, . . . , fn, encontrar p∈Pn, con p(xi) = fi,i= 0, . . . , n. Sabemos que este problema tiene una única solución. Para resolverlo expresamos p(x) = Pn i=0 cjxjy reducimos el problema al sistema Vc=f, donde V=V(x0, . . . , xn) := M1, x, . . . , xn x0, x1, . . . , xn=          1x0· · · xn 0 1x1· · · xn 1 . . .. . . 1xn· · · xn n          es la matriz de Vandermonde en los nodos x0, . . . , xn,c= (c0, c1, . . . , cn)T,f= (f0, f1, . . . , fn)T. Nuestro objetivo es estudiar la resolución del sistema Vc=f, es decir, determinar los coeficientes respecto de la base de monomios del interpolante. 6 La solución del problema de Lagrange puede expresarse a través de la fórmula de Lagrange p(x) = n X i=0 fjlj(x), lj(x) = Y k6=j x−xk xj−xk . Dado que la fórmula de Lagrange resuelve el problema de interpolación explícitamente en términos de los datos fdel segundo miembro, debe estar relacionado con la matriz inversa de V. Llamando l= (l0, . . . , ln)Ta la base de Lagrange y t= (t0, . . . , tn)Ta la base de monomios tj(x) := xj,j= 0, . . . , n, podemos comparar las expresiones del interpolante respecto a ambas bases lTf=tTc. Utilizando la relación c=V−1fobtenemos lTf=tTV−1f, y, como esta relación debe verificarse para todo f, deducimos que lT=tTV−1, es decir, la matriz de cambio de base de la base de Lagrange respecto de la base de monomios es la matriz inversa de Vandermonde. Considerando que Y k6=j (x−xk) = xn−X k6=j xixn−1+· · · + (−1)pX k1<···<kp∈{0,...,n}\{j} xk1· · · xkpxn−p+· · · + (−1)nY k6=j xk, tenemos que el coeficiente de ljen xies de la forma v(−1) ij =(−1)iP#K=n−i,K⊆{0,...,n}\{j}Qk∈Kxk Qk6=j(xj−xk),(3.1) obteniéndose el elemento (i, j)de la inversa de la matriz de Vandermonde. 4. Fórmula de Newton y factorización LU La fórmula de Newton del interpolante polinómico p(x) = n X j=0 [x0, . . . , xj]fωj(x) 7 permite expresar el interpolante en términos del vector de diferencias divididas d= (d0, . . . , dn), dj:= [x0, . . . , xj],j= 0, . . . , n, y la base de Newton ω= (ω0, . . . , ωn), con ωj(x)=(x−x0)· · · (x−xj−1). Cada elemento ωjde la base de Newton es un polinomio mónico de grado jcon la propiedad ωj(xi)=0, si j > i. Esto implica que la matriz de colocación en los nodos Mω0,...,ωn x0,x1,...,xnes una matriz triangular inferior L:=              1 0 0 · · · 0 1x1−x00· · · 0 1x2−x0(x2−x0)(x2−x1).... . . . . .. . ....0 1xn−x0(xn−x0)(xn−x1)· · · (xn−x0)· · · (xn−xn−1)              , cuyo elemento (i, j)es lij =ωj(xi) = Qj−1 k=0(xi−xk),j≤i. Evaluando en cada punto la expresión del polinomio de interpolación p=ωTd, obtenemos fi=ωT(xi)d,i= 0, . . . , n, obteniéndose el sistema Ld=f. Es decir, el vector de diferencias divididas es la solución del sistema triangular Ld=f. Podemos calcular los elementos de la matriz Lpor recurrencia, porque los elementos de la columna j-ésima de la matriz Lpueden obtenerse de la columna anterior de la siguiente manera lij =li,j−1(xi−xj−1),(4.1) partiendo de li0= 1,i= 0, . . . , n. Si aplicamos la fórmula de Newton a cada monomio obtenemos (1, . . . , xn)=(ω0(x), . . . , ωn(x))M1, . . . , xn [x0],[x0, x1],...,[x0, . . . , xn]. Teniendo en cuenta que [x0, . . . , xi]xj= 0 si i>j, se deduce que la matriz de cambio de base entre la base de Newton y la base de monomios es la matriz triangular superior U:=              1x0x2 0· · · xn 0 01[x0, x1]x2· · · [x0, x1]xn 0 0 1 .... . . . . .. . ....[x0, . . . , xn−1]xn 0 0 · · · 0 1              . 8 Parece ser que para calcular uij := [x0, . . . , xi]xjes necesario realizar varias diferencias y restas. Sin embargo, podemos obtener el valor uij en términos de los x0, . . . , xiutilizando una relación en la que solo aparecen sumas de potencias de los nodos. De acuerdo con la regla de Leibniz para diferencias divididas, tenemos [x0, . . . , xi]xj=xi[x0, . . . , xi]xj−1+ [x0, . . . , xi−1]xj−1, es decir uij =ui−1,j−1+xiui,j−1.(4.2) Usando esas relaciones puede calcularse la fija i-ésima partiendo de la fila (i−1)-ésima, teniendo en cuenta que uii = 1,i= 0, . . . , n, y uij = 0,j < i. Se deduce la siguiente fórmula para la fila i-ésima en términos de la fila (i−1)-ésima uij =ui−1,j−1+xiui−1,j−2+x2 iui−1,j−3+· · · +xj−i iui−1,i−1, o equivalentemente [x0, . . . , xi]xj= j−i X k=0 xk i[x0, . . . , xi−1]xj−1−k, j ≥i. Por inducción se demuestra que uij = [x0, . . . , xi]xj=X α0+···+αi=j−i xα0 0· · · xαi i. Si en la relación de cambio de bases tT=ωTU, tomamos matrices de colocación se deduce que M1, x, . . . , xn x0, x1, . . . , xn=Mω0, . . . , ωn x0, x1, . . . , xnU, lo que implica que V=LU, (4.3) es decir, las matrices LyUforman la factorización de Crout de la matriz de Vandermonde, donde la matriz triangular superior Utiene unos en la diagonal. La factorización LU se utiliza con frecuencia para resolver el sistema Vc=f. Llamando dal vector de diferencias divididas, tenemos Ld=f, Uc=d, y la solución del sistema con matriz de Vandermonde se reduce a la resolución consecutiva de dos sistemas triangulares con matrices LyU. Estos sistemas intermedios relacionan la solución con un vector intermedio d, el vector de las diferencias divididas y, por tanto, están relacionados directamente con la fórmula de Newton de interpolación. 9 La siguiente factorización es: VN=˜ L˜ U=          1000 1100 1210 1331                   1 0 0 0 0 1/3 1/9 1/27 002/9 2/9 0 0 0 2/9          , y sus normas son: ||˜ L||∞= 8,|| ˜ U||∞= 1, ||˜ L−1||∞= 8,|| ˜ U−1||∞= 9. En la siguientes tablas podemos ver los condicionamientos de estas matrices. κ∞(VN)κ∞(L)κ∞(U)κ∞(L)κ∞(U) 216 416 104 4 κ∞(VN)κ∞(˜ L)κ∞(˜ U)κ∞(˜ L)κ∞(˜ U) 216 576 64 9 Orden de Leja Vamos a hacer un análisis similar al realizado con la ordenación natural. Tomamos n= 3, los puntos equidistantes siguiendo la ordenación de Leja son: ˜x0= 1,˜x1= 0,˜x2=1 3,˜x3=2 3. La matriz de Vandermonde en estos puntos es VL=V(˜x0,˜x1,˜x2,˜x3) =          1 1 1 1 1 0 0 0 1 1/3 1/9 1/27 1 2/3 4/9 8/27          . Análogamente al apartado anterior, vamos a dar dos factorizaciones LU distintas. La primera de ellas: VL=LU =          1 0 0 0 1−1 0 0 1−2/3−2/9 0 1−1/3−2/9−2/27                   1 1 1 1 0 1 1 1 0014/3 0 0 0 1          , 16 ||L||∞= 2,||U||∞= 4, ||L−1||∞= 36,||U−1||∞= 2,3333. La segunda factorización es: VL=˜ L˜ U=          1 0 0 0 1 1 0 0 1 2/310 1 1/311                   1 1 1 1 0−1−1−1 0 0 −2/9−8/27 0 0 0 −2/9          , ||˜ L||∞= 3,3333,|| ˜ U||∞= 4, ||˜ L−1||∞= 2,6667,|| ˜ U−1||∞= 22,5. Las tablas mostradas a continuación recogen los condicionamientos de estas matrices: κ∞(VL)κ∞(L)κ∞(U)κ∞(L)κ∞(U) 216 672 72 9,3333 κ∞(VL)κ∞(˜ L)κ∞(˜ U)κ∞(˜ L)κ∞(˜ U) 216 800 8,8889 90 9.2. K= [−1,1] El intervalo que vamos a considerar en esta sección es [−1,1] yn= 3. En este intervalo tomaremos las siguientes ordenaciones: 1. La ordenación natural. 2. La ordenación de Leja. 3. Los puntos ordenados crecientemente según su distancia a cero. Si el intervalo es [−1,1], esto equivale a ordenar según su distancia al centro del intervalo. Por ello, vamos a llamarla ordenación central. De la misma manera que en los ejemplos en [0,1], consideraremos dos factorizaciones LU distintas. Orden natural Tomando el orden natural tenemos que los nodos son: x0=−1, x1=−1 3, x2=1 3, x3= 1. 17 La matriz de Vandermonde sobre estos nodos es: VN=V(x0, x1, x2, x3) =          1−1 1 −1 1−1/3 1/9−1/27 1 1/3 1/9 1/27 1 1 1 1          . Una de las factorizaciones es la siguiente: VN=LU =          1 0 0 0 1 2/3 0 0 1 4/3 8/9 0 128/3 16/9                   1−1 1 −1 0 1 −4/3 13/9 0 0 1 −1 0 0 0 1          . A continuación escribimos las normas de estas matrices y de sus inversas: ||L||∞= 7,4444,||U||∞= 4, ||L−1||∞= 4,5,||U−1||∞= 2,4444. La segunda factorización: VN=˜ L˜ U=          1000 1100 1210 1331                   1−1 1 −1 0 2/3−8/9 26/27 0 0 8/9−8/9 0 0 0 16/9          , y sus normas son: ||˜ L||∞= 8,|| ˜ U||∞= 4, ||˜ L−1||∞= 8,|| ˜ U−1||∞= 3,0625. En las siguientes tablas podemos ver los condicionamientos de estas matrices. κ∞(VN)κ∞(L)κ∞(U)κ∞(L)κ∞(U) 18 327,5556 33,5 9,7778 κ∞(VN)κ∞(˜ L)κ∞(˜ U)κ∞(˜ L)κ∞(˜ U) 18 784 64 12,25 18 Orden de Leja Siguiendo el mismo esquema que el utilizado con el orden natural vamos a ver el orden de Leja. Los puntos equidistantes siguiendo el orden de Leja son: ˜x0= 1,˜x1=−1,˜x2=−1 3,˜x3=1 3. La matriz de Vandermonde en estos puntos es VL=V(˜x0,˜x1,˜x2,˜x3) =          1 1 1 1 1−1 1 −1 1−1/3 1/9−1/27 1 1/3 1/9 1/27          . Análogamente al apartado anterior, vamos a dar dos factorizaciones LU distintas. La primera de ellas: VL=LU =          1 0 0 0 1−2 0 0 1−4/3−8/9 0 1−2/3−8/9−16/27                   1 1 1 1 0 1 0 1 001−1/3 0 0 0 1          , ||L||∞= 3,222,||U||∞= 4, ||L−1||∞= 4,5,||U−1||∞= 3,3333. La segunda factorización es: VL=˜ L˜ U=          1 0 0 0 1 1 0 0 1 2/310 1 1/311                   1 1 1 1 0−2 0 −2 0 0 −8/9 8/27 0 0 0 −16/27          , ||˜ L||∞= 3,3333,|| ˜ U||∞= 4, ||˜ L−1||∞= 2,6667,|| ˜ U−1||∞= 3,1875. Podemos ver los condicionamientos de estas matrices en las siguientes tablas: 19 κ∞(VL)κ∞(L)κ∞(U)κ∞(L)κ∞(U) 18 193,3333 14,5 13,3333 κ∞(VL)κ∞(˜ L)κ∞(˜ U)κ∞(˜ L)κ∞(˜ U) 18 113,3333 8,8889 12,75 Orden central Por último, tomemos los puntos ordenados crecientemente según su distancia a cero: ¯x0=1 3,¯x1=−1 3,¯x2=−1,¯x3= 1. La matriz de Vandermonde en estos puntos es VL=V(¯x0,¯x1,¯x2,¯x3) =          1 1/3 1/9 1/27 1−1/3 1/9−1/27 1−1 1 −1 1 1 1 1          . La primera factorización es la siguiente: VC=LU =          1 0 0 0 1−2/3 0 0 1−4/3 8/9 0 1 2/3 8/9 16/9                   1 1/3 1/9 1/27 0 1 0 1/9 0 0 1 −1 0 0 0 1          , ||L||∞= 4,3333,||U||∞= 2, ||L−1||∞= 4,5,||U−1||∞= 2. A continuación escribimos la segunda factorización: VC=˜ L˜ U=          1 0 0 0 1 1 0 0 1 2 1 0 1−111                   1 1/3 1/9 1/27 0−2/3 0 −2/27 008/9−8/9 0 0 0 16/9          , 20 ||˜ L||∞= 4,|| ˜ U||∞= 1,7778, ||˜ L−1||∞= 8,|| ˜ U−1||∞= 1,6875. Comparemos los condicionamientos de estas matrices: κ∞(VC)κ∞(L)κ∞(U)κ∞(L)κ∞(U) 18 78 19,54 κ∞(VC)κ∞(˜ L)κ∞(˜ U)κ∞(˜ L)κ∞(˜ U) 18 96 32 3 10. Experimentos numéricos: factores triangulares En esta sección vamos a comparar normas y condicionamientos de las matrices de las dos factorizaciones LU propuestas para distintas ordenaciones. Para ello vamos a considerar dos intervalos distintos. Los cálculos se han realizado en doble precisión con MATLAB. 10.1. Intervalo [0,1] En este caso K= [0,1]. Denotemos por LU la factorización en la que Utiene unos en la diagonal y por ˜ L˜ Ula factorización en la que ˜ Ltiene unos en la diagonal. Vamos a considerar puntos equidistantes en el intervalo [0,1] con dos ordenaciones: 1. Orden natural 2. Orden de Leja Comparemos primero las normas de Ly˜ Lcon las dos ordenaciones. ||L||∞ n Natural Leja 32,8889 2 43,2188 2 53,5104 2,0112 94,4583 2,0425 19 6,1522 2,0498 ||˜ L||∞ n Natural Leja 3 8 3,3333 4 16 4 5 32 4,65 9 512 7,5726 19 5,2429 ×10511,7892 21 Vemos que, tanto para ||L||∞como para ||˜ L||∞, es mejor la ordenación de Leja. Hemos podido probar la siguiente proposición sobre ||˜ L||∞y||˜ L−1||∞: Proposición 2. Sea Vla matriz de Vandermonde en los nodos x0, x1, . . . , xnequidistantes en [0,1] con la ordenación natural, y la factorización ˜ L˜ U. Se tiene ||˜ L||∞= 2n ||˜ L−1||∞= 2n Demostración. Sea L= (lij)0≤i,j≤nla matriz resultante de la factorización LU con Ucon unos en la diagonal. Ya hemos visto que lij =ωj(xi), por tanto: ˜ lij =ωj(xi) ωj(xj)=Qj−1 k=0(xi−xk) Qj−1 k=0(xj−xk) equidistantes =(i−j+ 1)(i−j+ 2). . . (i−0) j(j−1). . . 1=i j Por tanto, ||˜ L||∞= m´ax i=0,...,n i X j=0 |˜ lij|= m´ax i=0,...,n i X j=0 i j= m´ax i=0,...,n2i= 2n Denotemos por l(−1) ij y por ˜ l(−1) ij los elementos de L−1y˜ L−1respectivamente. l(−1) ij =1 ω0 i(xj)=⇒˜ l(−1) ij =ω0 i(xi) ω0 i(xj)=Qi−1 k=0(xi−xk) Qk=0,...,i k6=j (xj−xk) equidistantes = (−1)i−ji j Así, ||˜ L−1||∞= m´ax i=0,...,n i X j=0 |˜ l(−1) ij |= m´ax i=0,...,n i X j=0 i j= m´ax i=0,...,n2i= 2n De la proposición anterior se deduce el valor del condicionamiento de ˜ L. Por tanto, κ(˜ L) = 22n. Veamos qué ocurre para las matrices triangulares superiores. ||U||∞ n Natural Leja 324 42,55 53,2 6,5200 910,2469 36,2812 19 222,5312 2,1332 ×103 || ˜ U||∞ n Natural Leja 3 1 4 4 1 5 5 1 6 9 1 10 19 1 20 22 En contraste con lo que sucede con ||L||∞y||˜ L||∞, vemos que para las normas de ||U||∞ y|| ˜ U||∞es preferible la ordenación natural. Las siguientes tablas nos muestran los condicionamientos de estas matrices triangulares de las dos factorizciones. Natural Leja Natural Leja nκ∞(L)κ∞(U) 3 104 72 4 9,3333 4549,3333 341,3333 6,2500 16,875 52,9253 ×1031,6760 ×10311,5200 33,6432 92,4370 ×1061,1165 ×106138,6496 684,7808 19 5,2495 ×1013 1,7479 ×1013 1,0119 ×1051,5611 ×106 Natural Leja Natural Leja nκ∞(˜ L)κ∞(˜ U) 3 64 8,8889 9 90 4 256 16 26,6667 480 51,0240 ×10314,88 88,5417 4,0625 ×103 92,6214 ×10546,1569 1,3840 ×1049,0002 ×106 19 2,7488 ×1011 82,3227 7,1536 ×1097,6698 ×1015 Vemos que κ∞(L)es mejor con la ordenación de Leja aunque con el orden natural se obtienen resultados similares. En cambio, los mejores resultados para κ∞(U)se obtienen con el orden natural. Para la factorización ˜ L˜ Uobtenemos conclusiones análogas. Con el orden de Leja κ∞(˜ L)es más bajo, pero κ∞(˜ U)es mejor con el orden natural. El hecho de que κ∞(L)yκ∞(˜ L)sea menor para el orden de Leja es el esperado ya que es bien conocido que el pivotaje parcial controla el tamaño de los elementos de la matriz triangular inferior. Como ya hemos visto, el orden de Leja corresponde a seguir una estrategia de pivotaje parcial. En [9] se justifica el buen condicionamiento de ˜ U, y por tanto de U, para nodos positivos con la ordenación natural. En este artículo se prueba que si existe una estrategia de pivotaje óptima para reducir el condicionamiento de Skeel de la matriz triangular superior U, entonces esta estrategia coincide con el pivotaje parcial escalado para una norma ||.|| estrictamente monótona y en [6] se prueba que dicha estrategia aplicada a la eliminación gaussiana de una matriz TP no da lugar a cambio de filas. Una norma ||.|| se dice estrictamente monótona si, para cualesquiera vectores u= (u1, . . . , un),v= (v1, . . . , vn)con |uj|≥|vj|,∀j= 1, . . . , n 23 entonces ||u|| ≥ ||v|| y si además para algún j|uj|>|vj|entonces ||u|| >||v||. Así, los resultados mencionados dan una justificación teórica del buen condicionamiento de Uy˜ U para el orden natural. Para terminar esta sección, vamos a analizar los condicionamientos de las matrices de las dos factorizaciones LU con la ordenación natural y de Leja de los puntos de Chebyshev en [0,1]. La matriz de Vandemonde está mejor condicionada para puntos de Chebyshev como se muestra en el capítulo 21 de Higham [8] donde se recopilan resultados sobre los condicionamientos de matrices de Vandermonde con distintas distribuciones de nodos. Las dos siguientes tablas muestran los condicionamientos de las matrices de las dos factorizaciones. Natural Leja Natural Leja nκ∞(L)κ∞(U) 3112,5004 80,4374 4,1537 8,8590 4512,6342 327,1301 6,3730 20,0671 52,2857 ×1031,3061 ×10311,0868 26,1224 97,9909 ×1053,3996 ×105124,4100 642,6657 19 1,2943 ×1012 3,5654 ×1011 7,7281 ×1041,5873 ×106 Natural Leja Natural Leja nκ∞(˜ L)κ∞(˜ U) 353,4558 9,3137 10,7657 94,4402 4158,0263 17,7082 32,4371 568,3153 5588,4486 24,9545 106,5281 2,1546 ×103 98,0701 ×10447,3294 2,9425 ×1042,6893 ×106 19 2,2591 ×1010 149,7096 2,2093 ×1011 8,7018 ×1013 Se obtienen los mismos resultados que para puntos equidistantes. κ∞(L)es mejor con el orden de Leja mientras que κ∞(U)lo es con el orden natural. Ocurre lo mismo para la factorización ˜ L˜ U:κ∞(˜ L)es mucho más bajo con Leja y con el orden natural se obtienen los mejores resultados para κ∞(˜ U). 10.2. Intervalo [−1,1] Ahora el intervalo que vamos a considerar es el [−1,1]. Vamos a tomar los puntos equidistantes en este intervalo con diferentes ordenaciones: 1. La ordenación natural, es decir, los puntos ordenados de menor a mayor. 24 2. La ordenación de Leja. 3. La ordenación central, es decir, los puntos ordenados crecientemente según su distancia a cero. En este caso también vamos a comparar las dos factorizaciones: LU y˜ L˜ U. La primera se caracteriza por ser Ula que tiene unos en la diagonal y en la segunda es la matriz triangular inferior, ˜ L, la que tiene unos en la diagonal. Empecemos comparando las normas de las matrices triangulares inferiores: ||L||∞ n Natural Leja Central 37,4444 3,2222 4,3333 410,5 3,625 5,75 514,3408 3,8032 6,5392 92,4315 4,07726 8,3627 19 429,8239 4,1264 12,9618 ||˜ L||∞ n Natural Leja Central 3 8 3,3333 4 4 16 4 8 5 32 4,65 13 9 512 7,5762 87 19 5,2429 ×10511,7892 9,899 ×103 Con el orden de Leja es con el que se obtienen las mejores normas para ambas factorizaciones. Esto es debido, de nuevo, a que el orden de Leja es similar a hacer pivotaje parcial y el pivotaje parcial controla la norma de la matriz triangular inferior. Las siguientes tablas muestran las normas de Uy˜ U. ||U||∞ n Natural Leja Central 3 4 4 2 46,125 5 2 59,0416 62,24 942,9016 10 3,6214 19 2,3971 ×10349,5469 9,0455 || ˜ U||∞ n Natural Leja Central 3 4 4 1,7778 4 5 5 1,5 5 6 6 1,2499 9 10 10 1,125 19 20 20 1,0556 En el caso de las matrices Uy˜ Ula ordenación que da mejores resultados sobre las normas es la central. Además, vemos que la norma de ˜ Udisminuye cuando la dimensión del problema crece. Comparemos los condicionamientos de estas matrices con las diferentes ordenaciones: natural, Leja y central. Primero veamos los de la factorización LU en la tabla mostrada a continuación: 25 Una de las conclusiones que hemos obtenido es que, dada una matriz de Vandermonde V basada en los nodos x0, . . . , xntales que xi≥0,i= 0, . . . , n podemos realizar su factorización LU (asociada a la fórmula de Newton), de modo que el cálculo de L,L−1,UyU−1puede realizarse con alta precisión relativa. En los experimentos numéricos hemos tomado dos factorizaciones distintas de la matriz de Vandermonde V. Una de las factorizaciones tiene unos en la diagonal de Uy la hemos denotado por LU. La factorización LU está relacionada de forma natural con la fórmula de interpolación de Newton. La segunda factorización analizada, ˜ L˜ U, se caracteriza por ser ˜ Lla que tiene unos en la diagonal y está asociada a la eliminación gaussiana. En el primer caso estudiado hemos tomado los puntos en el intervalo [0,1] con dos ordenaciones distintas: natural y de Leja. Con la factorización LU, el producto de condicionamientos es menor para la ordenación natural, aunque κ∞(L)sea mejor con Leja. La factorización ˜ L˜ U con el orden de Leja da lugar al producto de condicionamientos ligeramente más bajo. Sin embargo, no existen diferencias significativas cuando utilizamos el orden natural con la factorización LU y ésta última es más natural dada su relación con la fórmula de Newton y con las diferencias divididas. En el caso del intervalo [−1,1], además del orden natural y el de Leja, también hemos propuesto otra ordenación llamada orden central. Hemos visto que para la factorización LU es mejor la ordenación central y para ˜ L˜ Use obtienen los mejores resultados con el orden de Leja. Los productos de los condicionamientos son casi idénticos en ambos casos por lo que escogemos la factorización LU de la ordenación central por la misma razón que en el caso [0,1]: por su relación con las diferencias divididas, con la fórmula de Newton y porque es adecuada para la evaluación del polinomio interpolante. 32 12. Referencias [1] L. Bos, S. de Marchi, A. Sommariva, M. Vianello:Computing multivariate Fekete and Leja points by numerical linear algebra, SIAM J. Numer. Anal., 48 (2010), pp. 19841999. [2] J. M. Carnicer, T. N. T. Goodman, J. M. Peña:Roundoff errors for polynomial evaluation by a family of formulae, Computing, 82 (2008), pp. 199-215. [3] C. W. Cryer:Some properties of totally positive matrices, Linear Algebra and its Applications, 15 (1976), pp. 1-25. [4] C. de Boor, A. Pinkus:Backward error analysis for totally positive linear systems, Numer. Math, 27 (1976/77), pp. 485-490. [5] J. Demmel, M. Gu, S. Eisenstat, I. Slapni˘ nar, K. Veselić, Z. Drma˘ c:Computing the singular value decomposition with high relative accuracy, Lineal Algebra and its Applications, 299 (1999), pp. 21-80. [6] M. Gasca, J. M. Peña:Scaled pivoting in Gauss and Neville elimination for totally positive systems, Applied Numerical Mathematics, 13 (1993), pp. 354-355. [7] N. J. Higham:Stability analysis of algorithms for solving confluent Vandermonde-like systems, SIAM J. Matrix Anal. Appl.,1 (1990), pp. 23-41. [8] N. J. Higham:Accuracy and Stability of Numerical Algorithms, SIAM, 1996. [9] J. M. Peña:Pivoting strategies leading to small bounds of the errors for certain linear systems, IMA Journal of Numerical Analysis, 16 (1996), pp. 141-153. [10] L. Reichel:Newton interpolation at Leja points, BIT 30, 2 (1990), pp. 332-346. 33