Source: https://es.scribd.com/doc/56810937/Apuntes-de-Analisis-Numerico
Timestamp: 2016-07-28 15:17:25+00:00

Document:
UploadSign inJoinBooksAudiobooksComicsSheet MusicEditors' Picks BooksHand-picked favorites from our editorsEditors' Picks AudiobooksHand-picked favorites from our editorsEditors' Picks ComicsHand-picked favorites from our editorsEditors' Picks Sheet MusicHand-picked favorites from our editorsTop BooksWhat's trending, bestsellers, award-winners & moreTop AudiobooksWhat's trending, bestsellers, award-winners & moreTop ComicsWhat's trending, bestsellers, award-winners & moreTop Sheet MusicWhat's trending, bestsellers, award-winners & moreCategoriesArts & IdeasBiography & MemoirBusiness & LeadershipChildren'sComputers & TechnologyCooking & FoodCrafts & HobbiesFantasyFiction & LiteratureHappiness & Self-HelpHealth & WellnessHistoryHome & GardenHumorLGBTMystery, Thriller & CrimePolitics & EconomyReferenceReligionRomanceScience & NatureScience FictionSociety & CultureSports & AdventureTravelYoung AdultCategoriesArts & IdeasBiography & MemoirBusiness & LeadershipChildren'sComputers & TechnologyCooking & FoodFantasyFiction & LiteratureHappiness & Self-HelpHealth & WellnessHistoryHome & GardenHumorLGBTMystery, Thriller & CrimePolitics & EconomyReferenceReligionRomanceScience & NatureScience FictionSociety & CultureSports & AdventureTravelYoung AdultCategoriesAdaptationsChildren’sCrime & MysteryFictionHumorMangaNonfictionRomanceSciFi, Fantasy & HorrorSuperheroesYoung AdultPublishersArcanaArchie ComicsBOOM! StudiosDynamiteIDW PublishingKingstone ComicsMarvel ComicsSpace Goat ProductionsTop Cow ComicsValiant ComicsZenescopeDifficultyBeginnerIntermediateAdvancedMixedInstrumentBrassDrums & PercussionGuitar, Bass, and FrettedPianoStringsVocalWoodwindsGenreClassicalCountryFolkJazz & BluesMovies & MusicalsPop & RockReligious & HolidayStandardsWelcome to Scribd! Start your free trial and access books, documents and more.Find out moreANÁLISIS NUMÉRICOMiguel Alemán Flores, Luis Álvarez León y Javier Sánchez Pérez Departamento de Informática y Sistemas Universidad de Las Palmas Campus de Taﬁra 35017 Las Palmas, España Email: {maleman,lalvarez,jsanchez}@dis.ulpgc.es
A continuación, mostraremos un resultado que indica el error de redondeo máximo que se produce al aproximar un número real cualquiera en una aritmética de precisión ﬁnita. Teorema 2 Sean ymin , ymax los valores positivos menor e e y mayor de una aritmética de precisión ﬁnita. Sea u la unidad de redondeo de dicha aritmética. Si un número real z veriﬁca que ymin <| z |< ymax , entonces e e donde z es el número más cercano a z en la aritmética. e | z − z |≤| z | u e
que las variables A y B están cercanas entre sí con una tolerancia T OL si se cumple que | A − B |≤ max {| A |, | B |} T OL Este criterio es simétrico en el sentido de que trata de igual modo los números A y B. También se puede utilizar un criterio más simple, como | A − B |≤| A | T OL pero en este caso le estamos dando una signiﬁcación especial a A con respecto a B. Estos criterios de comparación de números funcionan bien salvo cuando los números A y B están muy próximos a 0. Por ejemplo, si B = 0, los criterios anteriores quedan | A |≤| A | T OL lo cual es imposible (si T OL < 1), salvo que A también sea 0. Para evitar este comportamiento, se puede añadir al criterio un valor > 0 de la siguiente forma: | A − B |≤ (max {| A |, | B |} + ) T OL
donde a0 = 1 y, en general, an = 0 o an = 1. Además, para un número natural t cualquiera tenemos que Ã t ! t X an X an 1 e e 2 ≤z≤2 + t 2n 2n 2 n=1 n=1
Demostración. Un número real cualquiera z, que tomaremos positivo sin pérdida de generalidad, se puede expresar como ∞ X an z = 2e 2n n=1
Por el problema anterior, el número que está a la derecha de la desigualdad también pertenece a la aritmética de precisión ﬁnita, y por tanto | z − z |≤ e 2e 2−t 2
Programa 4 Programa en Fortran 77 que determina si dos variables A, B son iguales con una tolerancia T OL (tomando el máximo de A, B), con = 10−10 READ *,A,B,TOL IF(IGUAL(A,B,TOL).EQ.0) THEN PRINT *,’A=B segun la tolerancia TOL’ STOP ELSE PRINT *,’A distinto de B segun la tolerancia TOL’ STOP ENDIF END FUNCTION IGUAL(A,B,TOL) IF(ABS(A).GT.ABS(B)) THEN IF(ABS(A-B).LE.(TOL*(ABS(A)+10.**(10.))) THEN IGUAL=0 RETURN ELSE IGUAL=1 RETURN ENDIF ELSE IF(ABS(A-B).LE.(TOL*(ABS(B)+10.**(10.))) THEN IGUAL=0 RETURN ELSE IGUAL=1
Ahora bien, como a0 = 1, se tiene que 2e < 2 | z | y, por tanto | z − z |≤| z | 2−t =| z | u e con lo que queda demostrado el teorema. Problema 9 (2 puntos) Dado un número z = e P an 2e t n=1 2n , en una aritmética de precisión ﬁnita. Calcular el número inmediatamente inferior a él en dicha aritmética. Un resultado importante para la comparación de dos números es el siguiente: Teorema 3 Si z1 , z2 ∈ A son distintos entonces e e Demostración: Ejercicio En muchos algoritmos, el test de parada incluye el hecho de que dos variables estén próximas entre sí. para ello se ﬁja un umbral o tolerancia T OL que por supuesto será mayor que la unidad de redondeo u y expresaremos | z1 − z2 |≥ max {| z1 |, | z2 |} u e e e e
RETURN ENDIF ENDIF END 1 Asociado a cualquier aritmética de precisión ﬁnita de números reales, existen 4 operaciones básicas, que son la suma, la resta, la multiplicación y la división de números reales dentro de la aritmética. Nosotros no vamos a entrar en este curso en cómo se pueden deﬁnir algorítmicamente estas operaciones. Solamente queremos mencionar que, a menudo, para minimizar el efecto de los redondeos en las operaciones, antes de realizarlas se aumenta la precisión de los números reales (por ejemplo pasando de simple precisión a doble precisión) para, a continuación, realizar la operación en una aritmética de mayor precisión y, ﬁnalmente, el resultado se redondea para pasarlo a la precisión inicial.
C=0.01 D=0 DO 1 K=1,2**7 B=B+A CONTINUE DO 2 K=1,100 D=D+C CONTINUE PRINT *,(1-B)*(10**10) PRINT *,(1-D)*(10**10) END
Además, este programa permite identiﬁcar la base de la aritmética con la que trabaja el ordenador. Como conclusión de este apartado, podemos extraer que, para ser más precisos numéricamente, cuando trabajamos con números más pequeños que la unidad deberíamos pensar en términos de 2−m en lugar de 10−m , que es como solemos hacerlo. Errores por Cancelación. Estos errores se producen al restar números de aproximadamente la misma magnitud. Hay que tener en cuenta que, al realizar operaciones sobre una variable, los errores de redondeo se van acumulando en la parte menos signiﬁcativa del número (los dígitos de menos valor), dejando relativamente intacta la parte más signiﬁcativa del número, que corresponde a los dígitos de mayor valor. Por ello, al restar dos números de magnitud parecida, se cancelan las partes signiﬁcativas, quedando la aportación de los dígitos de menos valor, que es donde más error hay. Por ejemplo, en el programa Fortran anterior, se ha utilizado este fenómeno de cancelación para poner de maniﬁesto la diferencia entre trabajar con bases distintas. En los algoritmos, muchas veces se intenta evitar la posibilidad de restar 2 números que pudieran ser de magnitud parecida. Por ejemplo, en la conocida fórmula del cálculo de raíces de un polinomio de grado 2, ax2 + bx + c = 0 (con a 6= 0) √ −b ± b2 − 4ac x= 2a una forma de evitar la cancelación que se produce cuando √ b ≈ b2 − 4ac consiste en calcular primero la raíz de mayor valor absoluto, es decir √ ¢ ¡ − b + sign(b) b2 − 4ac x1 = 2a y después la segunda raíz x2 utilizando la relación x1 x2 = c a. Por lo tanto, en los algoritmos, se deberá evitar, en la medida de lo posible, la resta de variables que tengan una magnitud cercana. Problema 10 (1 punto) Calcular las raíces del polinomio P (x) = x2 − 2x + 0.01 evitando los errores de cancelación.
Fuentes de errores numéricos Dentro de las posibles fuentes de errores numéricos, destacaremos 3 tipos: Errores de redondeo. Son los que se producen al ”redondear” un número real para poder expresarlo en una aritmética de precisión ﬁnita. Como vimos en la sección anterior, este error está controlado por la denominada unidad de redondeo, u = 2−t , de tal forma que, al tomar un número real z y aproximarlo en la aritmética por el valor z ∈ A más próximo, el error de redondeo tiene la e expresión: | z − z |≤| z | u e Errores de cambio de base. Este tipo de errores se produce al realizar un cambio de base para representar un número real. Como vimos en la sección anterior, las aritméticas estándares de ordenador trabajan en base 2. Sin embargo, los humanos pensamos y razonamos en términos de números en base 10. Por ejemplo, números tan naturales para nosotros como 0.1 no pueden representarse de forma exacta en una aritmética en base 2. Esto quiere decir que, al representar 0.1 el ordenador, va a producir un pequeño redondeo, y este pequeño error de redondeo se puede ir propagando hasta producir errores apreciables. Por ejemplo, parece razonable pensar que, cuando sumamos 100 veces el número 0.01, el resultado sea exactamente 1, pero, no es así. Sin embargo, si sumamos 128 = 27 veces el número 2−7 , el resultado sí es exactamente 1. Este resultado se pone de maniﬁesto en el siguiente programa Fortran: Programa 5 Programa en Fortran 77 para comprobar la diferencia entre trabajar en base 10 y trabajar en base 2. A=2**(-7.) B=0
Problema 11 (3 puntos) Escribir el código en Fortran 77 para implementar el cálculo de las raíces de ax2 + bx + c = 0 evitando los errores de cancelación y teniendo en cuenta las diferentes opciones que aparecen cuando a 6= 0 y a = 0.
Método de Newton-Raphson Éste es, sin duda, uno de los métodos más importantes y útiles para el cálculo de raíces. Dada una aproximación inicial de la raíz x0 , se busca, a partir de x0 , una aproximación mejor x1 de la raíz, de la siguiente forma: Se sustituye la función f (x) por el valor de su desarrollo de Taylor centrado en x0 hasta el orden 1, es decir f (x) ≈ f (x0 ) + f 0 (x0 )(x − x0 ) que corresponde a un polinomio de grado 1, y a continuación se calcula x1 como el cero de este polinomio, es decir: f (x0 ) x1 = x0 − 0 f (x0 ) y por tanto, de forma general, obtenemos, a partir de x0 una secuencia xn de valores que van aproximando la raíz, deﬁnidos por f (xn ) xn+1 = xn − 0 f (xn ) A continuación veremos una aplicación de este método para calcular la raíz cuadrada √ un número positivo A, de teniendo en cuenta que si x = A, entonces f (x) = x2 − A = 0. Programa 6 Programa en Fortran 77 para calcular una aproximación de la raíz cuadrada de un número positivo A con una tolerancia T OL, y un número máximo de iteraciones N max . READ *,A,TOL,Nmax IF(A.LE.0) THEN PRINT *,’El numero A no es positivo’ STOP ENDIF X0=(1+A)/2. DO 1 K=1,Nmax X1=X0-(X0*X0-A)/(2.*X0) IF(IGUAL(X0,X1,TOL).EQ.0) THEN PRINT *,’LA RAIZ DE A ES’,X0 STOP ELSE X0=X1 ENDIF CONTINUE PRINT *,’No máximo de iterac. excedido’ END
CÁLCULO DE LOS CEROS DE UNA FUNCIÓN En esta sección vamos a estudiar algunos métodos para calcular los ceros de una función de una variable, f (x), esto es, los valores de x para los cuales f (x) = 0.
Método de la bisección Se considera un intervalo [a, b] donde la función f (x) cambia de signo, es decir f (a)·f (b) < 0. El método consiste en ir dividiendo el intervalo [a, b] por la mitad de la siguiente forma: Se toma el punto medio a+b . Si f ( a+b ) = 0 ya 2 2 hemos encontrado la raíz x = a+b . En caso contrario, 2 si f ( a+b ) · f (b) < 0 entonces hacemos a = a+b y volvemos 2 2 a subdividir el nuevo intervalo [a, b]. Si, por el contrario, f (a) · f ( a+b ) < 0, entonces hacemos b = a+b y volvemos 2 2 a empezar. Las sucesivas subdivisiones del intervalo [a, b] van aproximando la raíz.
Problema 12 (2 puntos) Calcular 2 iteraciones del algoritmo de la bisección para buscar un cero de la función f (x) = x2 − 2 en el intervalo [−2, 0]. Problema 13 (3 puntos) Escribir el código en Fortran 77 para implementar el método de la bisección
Método de la Regula-falsi (regla de lo falso) Este método es una variación del anterior en el sentido siguiente: En lugar de tomar el punto medio a+b del in2 tervalo, se considera el punto de intersección de la recta que pasa por los puntos (a, f (a)) y (b, f (b)) con el eje x. Es decir, en el razonamiento anterior, se sustituye el valor xm = a+b por el valor 2 xm b−a =a− f (a) f (b) − f (a)
Problema 14 (2 puntos) Calcular 2 iteraciones del algoritmo de la regula-falsi para buscar un cero de la función f (x) = x2 − 2 en el intervalo [0, 2]. Problema 15 (3 puntos) Escribir el código en Fortran 77 para implementar el método de la Regula-falsi.
de una forma analítica. En este caso, se sustituye el valor f 0 (xn ) en el algoritmo, por el valor f (xn ) − f (xn−1 ) xn − xn−1 que corresponde a una aproximación de f 0 (xn ). Para iniciar el algoritmo, son necesarias dos aproximaciones iniciales, x0 y x1. Problema 16 (1 punto) Calcular una iteración del método de Newton-Raphson para calcular un cero de la función f (x) = x3 − 3 partiendo de x0 = 1. Problema 17 (1 punto) Calcular una iteración del método de la secante para calcular un cero de la función f (x) = x3 − 3 partiendo de x0 = 0, x1 = 1. Problema 18 (3 puntos) Escribir un programa en Fortran 77 que implemente el método de la Secante utilizando reales de doble precisión. Los datos de entrada son las aproximaciones iniciales, x0 y x1 , el número máximo de iteraciones, N max, y la tolerancia, T OL, para determinar la igualdad de dos números.
Como veremos posteriormente, la elección de las fórmulas anteriores equivale a aproximar f (x) por la parábola que pasa por los puntos (xn−3 , f (xn−3 )) , (xn−2 , f (xn−2 )) y (xn−1 , f (xn−1 )), y calcular posteriormente las derivadas de dicha parábola. Programa 7 Programa en Fortran 77 donde se muestra un ejemplo de manejo de números complejos. IMPLICIT COMPLEX (C) CX=(-1,0) CY=CF(CX) PRINT *,CY END FUNCTION CF(CX) IMPLICIT COMPLEX(C) CF=SQRT(CX) END
Práctica 2 (Método de Müller, 4 horas) Implementar el método de Müller. Crear un programa en Fortran 77 que tenga como datos de entrada: las tres primeras aproximaciones de la raíz, x0 , x1 y x2 , el número máximo de iteraciones N max, y la tolerancia T OL, para determinar la igualdad entre dos números. La función a la que se le calculan los ceros se deﬁne en el propio cuerpo del programa. Utilizar el método para calcular los posibles ceros de las siguientes funciones: 1. f (x) = x2 + 1 2. f (x) = (x2 + 1)x 3. f (x) = ex − 1 4. f (x) = x − 2 5. f (x) = 1
Método de Müller Este método es de utilidad para calcular raíces complejas de funciones, como por ejemplo polinomios. Es una generalización del método de Newton-Raphson, en el sentido de que, en lugar de quedarnos con la parte lineal del desarrollo de Taylor de la función, nos quedamos con los términos hasta el orden 2, de tal forma que hacemos f 00 (xn−1 ) f (x) ≈ f (xn−1 )+f (xn−1 )(x−xn−1 )+ (x−xn−1 )2 2
donde xn−1 es una aproximación de una raíz compleja de la función f (x). Para obtener una aproximación xn mejor de la raíz calculamos los ceros del polinomio de segundo grado anterior, es decir q −f 0 (xn−1 ) ± (f 0 (xn−1 ))2 − 2f (xn−1 )f 00 (xn−1 ) xn = xn−1 + f 00 (xn−1 ) De las dos posibles raíces, nos quedamos con aquélla que sea más cercana a xn−1 . Dicha raíz será la aproximación xn de la raíz de f (x) en la etapa n. En el caso en que f 00 (xn−1 ) = 0, calculamos xn por el método de Newton-Raphson. En el caso en que no conozcamos analíticamente el valor de la primera y segunda derivada de f (x), podemos utilizar las siguientes aproximaciones: f 00 (xn−1 ) ≈ 2
Nota: Utilizar como tolerancia T OL = 0.0001 y N max = 100. Para el ejemplo 1 tomar como datos iniciales x0 = (3, 0), x1 = (2, 0) x2 = (1, 0). Para el ejemplo 2 tomar como datos iniciales x0 = (3, 0), x1 = (2, 0), x2 = (1, 0) y x0 = (1, 0), x1 = (0.1, 0), x2 = (0.01, 0). Para el ejemplo 3 tomar como datos iniciales x0 = (3, 0), x1 = (2, 0), x2 = (1, 0). Para el ejemplo 4 tomar como datos iniciales x0 = (3, 0), x1 = (2, 0), x2 = (1, 0). Para el ejemplo 5 tomar como datos iniciales x0 = (3, 0), x1 = (2, 0), x2 = (1, 0).
Cálculo de las raíces de un polinomio Los polinomios son un tipo particular de funciones que, por su gran utilidad, requieren un análisis algo más detallado. Nos ocuparemos sólo de las raíces reales de los polinomios, aunque también hay que indicar que existen algoritmos versátiles para el cálculo de las raíces complejas, como, por ejemplo, el método de Müller, visto anteriormente. A menudo, los alumnos pueden tener la impresión de que los algoritmos y técnicas que se aprenden en una asignatura como análisis numérico les serán de poca utilidad en el futuro. Mi experiencia como docente en esta disciplina es que, con frecuencia, una vez terminada la carrera y en el desarrollo de la actividad profesional, aparecen problemas que, para su resolución, requieren el uso de alguna de las técnicas presentadas en esta asignatura. El siguiente ejemplo es una buena prueba de ello. Ejemplo 2 Actualmente están muy de moda los planes de pensiones. Las entidades ﬁnancieras venden a sus clientes los planes de pensiones de la siguiente forma, por ejemplo: si usted aporta durante 30 años 100.000 pesetas todos los años, aportación que se va incrementando cada año en un 10%, es decir el primer año 100.000, el segundo año 110.000, etc., entonces, le aseguramos que al ﬁnal del trigésimo año tendrá a su disposición la cantidad de 26.000.000 de pesetas. Ahora bien, el dato más importante para el futuro pensionista (que a menudo oculta la entidad ﬁnanciera) es el interés nominal anual que se está aplicando año tras año al dinero depositado. Si llamamos i al interés nominal anual que se aplica al dinero, la ecuación que debemos resolver para obtener i es
resultado muestra una forma rápida y sencilla de evaluar simultáneamente un polinomio y su derivada. Teorema 4 (Método de Horner). Sea P (x) = an xn + an−1 xn−1 + ...... + a0 , si deﬁnimos bk como bn = an bk = ak + bk+1 x0 entonces se veriﬁca que P (x0 ) = b0 P 0 (x0 ) = bn xn−1 + bn−1 xn−2 + ........ + b1 0 0 Demostración Sea el polinomio Q(x) = bn xn−1 + bn−1 xn−2 + ..... + b1 . Veamos que se veriﬁca que P (x) = (x − x0 )Q(x) + b0 Efectivamente, dado que ak = bk − bk+1 x0 y an = bn , obtenemos la igualdad anterior teniendo en cuenta que (x − x0 )Q(x) + b0 = bn xn + (bn−1 − bn x0 )xn + ..... + (b0 − b1 x0 ) Por último, obtenemos P 0 (x) = (x − x0 )Q0 (x) + Q(x) de donde sale obviamente que P 0 (x0 ) = Q(x0 ). Este teorema permite calcular el polinomio y su derivada en un punto de forma muy sencilla, como muestra el siguiente programa Fortran. Programa 8 El siguiente programa en Fortran 77 calcula la evaluación de un polinomio y su derivada en un punto X, almacenándolos en las variables P X y P P X. PARAMETER(NMAX=1000) DIMENSION A(0:NMAX) COMMON/POL/PX,PPX PRINT *,’Escribir Grado del Polinomio’ READ *,N IF(N.GT.NMAX) THEN PRINT *, ’Grado Superior al Maximo’ STOP ENDIF PRINT *,’EscriBir Coef. Polin.’ DO 1 K=0,N READ *,A(K) PRINT *,’Escribir valor de X’ READ *,X CALL HORNER(N,A,X) PRINT *,’P(X)= ’,PX PRINT *,’P ‘(X)= ’,PPX END
(100.000) (1.1) (1. + i)30−n = 26.000.000
Ahora bien, para calcular i, debemos calcular las raíces del polinomio en i dado por P (i) =
(100.000) (1.1) (1. + i)30−n − 26.000.000
El cálculo de las raíces de este polinomio nos lleva a i = 4.487%. Este ejemplo muestra como un problema ﬁnanciero sencillo nos lleva a la necesidad de calcular los ceros de un polinomio.
Algoritmo de Horner para evaluar un polinomio en un punto Dado un polinomio P (x) = an xn + an−1 xn−1 + ...... + a0 , éste se puede expresar también de la forma siguiente: P (x) = a0 + x (a1 + x (a2 + x (a3 + x(..... + x (an−1 + xan ))))) . Además, si queremos utilizar un método de cálculo de raíces como el de Newton-Raphson, necesitamos evaluar tanto el polinomio como su derivada. El siguiente
SUBROUTINE HORNER(N,A,X) DIMENSION A(0:*) COMMON/POL/PX,PPX PX=A(N) PPX=A(N) DO 1 K=N-1,1,-1 PX=PX*X+A(K) PPX=PPX*X+PX CONTINUE PX=PX*X+A(0) END
Demostración [Is-Ke] Pg. 126. Para la estimación del número de raíces reales negativas, se aplica el teorema anterior cambiando x por −x. Ejemplo 3 Sea P (x) = 3x4 + 10x3 − 10x − 3, los signos de los coeﬁcientes son: + + −−. Por tanto, hay un único cambio de signo y hay una raíz positiva. Si cambiamos x por −x, los signos de los coeﬁcientes son + − +−. Por tanto, hay 3 cambios de signo y hay una o tres raíces negativas. En este caso, las raíces son x = 1, −1, −3, − 1 . 3 Problema 19 (1 punto) Calcular una iteración del método de Müller para calcular un cero de la función f (x) = x3 − 3 partiendo de x0 = 1 (Calculando las derivadas de la función de forma exacta) y quedándonos con la raíz más cercana a x0 .
Nota: La declaración P ARAM ET ER(N M AX = 1000) permite deﬁnir constantes. La declaración DIM EN SION A(0 : N M AX) deﬁne un vector A, de reales en precisión simple, de tamaño N M AX + 1, y numerados desde 0 hasta N M AX. La declaración COM M ON/P OL/P X, P P X deﬁne la zona de memoria denomina P OL donde se encuentran las variables globales P X, P P X. Para que una subrutina pueda hacer uso de esas variables, debe incluir en su inicio la misma sentencia COM M ON. Otros resultados interesantes de utilidad para localizar en qué zonas pueden estar las raíces del polinomio son: Teorema 5 Sea un polinomio P (x) = an xn +an−1 xn−1 + ...... + a0 con an 6= 0, entonces las raíces reales de P (x) están en el intervalo ∙ ¸ maxk=0,..,n−1 | ak | maxk=0,..,n−1 | ak | −1 − ,1 + | an | | an | Demostración Veamos que si |x| > 1 + entonces |P (x)| > 0. Efectivamente, |P (x)| ≥ |an xn | −
k=0,,,n−1 maxk=0,..,n−1 |ak | , |an |
Problema 20 (2 puntos) Dado el polinomio P (x) = 2x3 + 3x2 + 4x + 5, evaluar el polinomio y su derivada en el punto x = 2, utilizando el algoritmo de Horner.
Problema 21 (1 punto) Calcular el número máximo de raíces positivas y negativas del polinomio x5 −35x3 +30x2 + 124x − 120, y localizarlas en un intervalo. Teorema 7 Entre dos raíces de una función derivable f (x) hay una raíz de f 0 (x). Demostración Teorema de Rolle. Teorema 8 La derivada k − esima P k) (x) del polinomio ´ P (x) = an xn + an−1 xn−1 + ......a0 es P k) (x) = k! an n! n−k an−1 (n − 1) ! n−k−1 + +...+ak x x (n − k)! (n − k − 1)! 1
1 − |x| ≥ k=0,,,n−1 1 − |x| n |x| n = ≥ |an | |x| − max |ak | k=0,,,n−1 |x| − 1 |x|n (|an | (|x| − 1) − maxk=0,,,n−1 |ak |) = >0 |x| − 1 = |an | |x|n − max |ak |
Demostración Es inmediato, derivando sucesivamente el polinomio P (x). Los dos resultados anteriores permiten aislar las posibles raíces de P (x) de la forma siguiente: Si llamamos maxk=0,..,n−1 |ak | Pmax a 1 + , entonces las m raíces distintas |an | x1 < x2 < .... < xm de P (x) están intercaladas con las 0 raíces x0 < x0 < .... < x0 1 2 m−1 de P (x), es decir −Pmax ≤ x1 ≤ x0 ≤ x2 ≤ x0 ≤ ... ≤ x0 1 2 m−1 ≤ xm ≤ Pmax Volviendo a aplicar este razonamiento sucesivamente sobre P 0 (x), P 00 (x), etc., para intercalar los ceros de una derivada con los ceros de la siguiente, podemos deducir el siguiente algoritmo para aislar todas las raíces de un Polinomio P (x):
Teorema 6 Sea un polinomio P (x) = an xn +an−1 xn−1 + ...... + a0 , entonces el número de raíces positivas es igual al número de cambios de signo en los coeﬁcientes an , ......, a0 (saltando los posibles coeﬁcientes nulos), o bien ese mismo número menos un número par.
1. Se parte del intervalo [−Pmax , Pmax ]
2. Se calcula la raíz x1 del Polinomio P n−1) (x) (que es un polinomio de grado 1) 3. Para k = n − 2, ..., 1 −Pmax <
-2.5 -1.25
50 25 0 0 -25 -50 -75 -100 1.25 2.5 x
Al ﬁnal del procedimiento, habremos aislado completamente a las raíces de P (x). Este procedimiento se puede utilizar para grados relativamente pequeños (n < 30), puesto que su utilización requiere el cálculo de factoriales, que se dispara rápidamente. Por ejemplo, 30! = 2. 6 × 1032 . Existen métodos mejores para el cálculo de raíces de polinomios, pero que utilizan técnicas más complejas. El método presente en el siguiente programa, que combina el aislamiento de las raíces del polinomio a través de los ceros de sus derivadas con el método de Newton-Raphson, funciona razonablemente bien para grados de polinomios pequeños. En el caso de raíces múltiples los resultados acumulan mayores errores de redondeo debido a que tanto el polinomio como su derivada son cero en el mismo punto. Ejemplo 4 Consideremos el polinomio P (x) = x4 − x3 − 7x2 +x+6, que tiene por raices x = 1, 3, −1, −2. Para este polinomio, tenemos que Pmax = 8. Por tanto, las raíces están en el intervalo [−8, 8]. Por otro lado su gráﬁca es
Polinomio P 0 (x) = 4x3 − 3x2 − 14x + 1 La derivada segunda de este polinomio es P 00 (x) = 12x2 − 6x − 14,cuyas raíces son x = −0.858, 1. 358 y cuya gráﬁca es
25 0 -2.5 -1.25 0 1.25 2.5 x
La derivada tercera de este polinomio es P 000 (x) = 24x − 6, cuya raíz es x = 0.25, y cuya gráﬁca es
20 50 0 -2.5 -1.25 0 1.25 2.5 x -2.5 -1.25 -25 25 0 0 1.25 2.5 x
La derivada de este polinomio es P 0 (x) = 4x3 − 3x2 − 14x + 1,cuyas raíces son x = −1. 574, 7. 05 × 10−2 , 2. 253 y cuya gráﬁca es
Polinomio P 000 (x) = 24x − 6 El método funcionaría de la siguiente forma: Primero calculamos el cero de P 000 (x), es decir x = 0.25, por tanto los ceros de P 00 (x) estarían en los intervalos [−2.166, 0.25] y [0.25, 2.166]. Puesto que hay cambio de signo de P 00 (x) en cada uno de estos intervalos, buscamos las raíces de P 00 (x) en esos intervalos, utilizando cualquier método numérico de los vistos anterioremente, por ejemplo, el método de la Regula-falsi, obteniendo −0.858 para el inter12
valo [−2.166, 0.25] y 1. 358 para el intervalo [0.25, 2.166]. Por tanto, las posibles raíces de P 0 (x) estarán en los intervalos [−4.5, −0.858], [−0.858, 1.358] y [1.358, 4.5]. Buscamos ahora las raíces de P 0 (x) es esos intervalos, obteniendo x = −1. 574, 7. 05 × 10−2 y 2. 253. Por tanto, los posibles ceros de P (x) estarán en los intervalos [−8, −1. 574], [−1.574, 7. 05×10−2 ], [7. 05×10−2 , 2. 253] y [2.253, 8]. Buscamos, ﬁnalmente, las raíces de P (x) en cada un de esos intervalos y obtenemos x = −2, −1, 1, 3. Problema 22 (2 puntos) Aislar en intervalos las raíces del polinomio P (x) = 20x3 − 45x2 + 30x − 1.
Programa 9 Programa en Fortran 77 donde se implementa la función ICEROP OL(A, R, T OL, N, N maxx), que devuelve las raíces reales de un polinomio. Dicha subrutina tiene como parámetros un vector A(), donde están los coeﬁcientes del polinomio, un vector R(), donde se guardan las raíces del polinomio una vez calculadas, la tolerancia T OL, con la que consideramos que dos números son iguales, el grado del polinomio N, y el número máximo de iteraciones N max xx, para el proceso de Newton-Raphson. También se deﬁne la función auxiliar RP (N, A, X1, X2, T OL, N maxx, R, L), que devuelve la raíz del polinomio que se obtiene aplicando el método de Newton-Raphson, tomando como valor inicial el punto medio del intervalo [X1, X2].
PARAMETER(NMAX=30) DIMENSION A(0:NMAX),R(0:NMAX-1) COMMON/POL/PX,PPX PRINT *,’Escribir Grado del Polinomio’ READ *,N IF(N.GT.NMAX) THEN PRINT *, ’Grado Superior al Maximo’ STOP ENDIF PRINT *,’Escribir Coef. Polin.’ DO 1 K=0,N 1 READ *,A(K) PRINT *, ’Escribir Tolerancia’ READ *,TOL PRINT *, ’Escribir No. iter. Max. para NewtonRaphson’ READ *,Nmaxx M=ICEROPOL(A,R,TOL,N,Nmaxx) PRINT *,’El Pol. tiene’,M,’ raices’ DO 9 K=0,M-1 9 PRINT *,R(K) END
DIMENSION A(0:*), R(0:*), F(0:NMAX), AP(0:NMAX), PI(0:NMAX+1) COMMON/POL/PX,PPX **** Calculo de los factoriales F(0)=1. DO 2 K=1,N 2 F(K)=F(K-1)*K *** Calculo intervalo inicial PMAX=ABS(A(0)) DO 3 K=1,N-1 IF(PMAX.LT.ABS(A(K)) THEN PMAX=ABS(A(K) ENDIF 3 CONTINUE PMAX=PMAX/ABS(A(N))+1. PI(0)=-PMAX PI(1)=-(A(N-1)*F(N-1))/(A(N)*F(N)) DO 10 K=2,N 10 PI(2)=PMAX *** Calculo de los coeﬁcientes del *** polinomio derivada DO 7 K=2,N 7 PI(K)=PMAX DO 4 K=N-2,0,-1 DO 5 L=0,N-K AP(L)=A(L+K)*(F(K+L)/F(L)) 5 CONTINUE ***CALCULAR LOS CEROS DE AP EN LOS INTERVALOS PI() DO 6 L=1,N-K PI(L)=RP(N-K,AP,PI(L1),PI(L),TOL,Nmaxx,R,L-1) 6 CONTINUE 4 CONTINUE *** Pasamos las raices al vector R() M=0 DO 8 K=1,N IF(R(K-1).EQ.0) THEN R(M)=PI(K) M=M+1 ENDIF 8 CONTINUE ICEROPOL=M END
FUNCTION ICEROPOL(A,R,TOL,N,Nmaxx) PARAMETER(NMAX=30)
FUNCTION RP(N,A,X1,X2,TOL,Nmaxx,R,L) DIMENSION A(0:*),R(0:*) COMMON/POL/PX,PPX R(L)=1. IF (X1.EQ.X2) THEN RP=X1 RETURN ENDIF RP=(X1+X2)/2. DO 1 K=1,Nmaxx CALL HORNER(N,A,RP) IF (PPX.EQ.0.) THEN IF(PX.EQ.0.) THEN
R(L)=0. RETURN ELSE RETURN ENDIF ELSE RP1=RP-PX/PPX IF(IGUAL(RP1,RP,TOL).EQ.0) THEN RP=RP1 R(L)=0. RETURN ELSE RP=RP1 ENDIF ENDIF CONTINUE END
Ejemplo 5 Consideremos una función f (x) = ex , vamos a interpolarla en los puntos x0 = 0, x1 = −1 y x2 = 1. Para calcular P2 (x), el polinomio interpolador de Lagrange en estos puntos, calcularíamos los polinomios base: (x + 1)(x − 1) −1 x(x − 1) 2 x(x + 1) 2
P2 (x) = e0 Problema 23 (2 puntos) Aislar en intervalos las raíces del polinomio P (x) = 2x3 + 3x2 − 12x + 1.
INTERPOLACIÓN DE FUNCIONES I El problema general de la interpolación de funciones consiste en, a partir del conocimiento del valor de una función (y eventualmente de sus derivadas) en un conjunto ﬁnito de puntos, aproximar el valor de la función fuera de ese conjunto ﬁnito de puntos.
Interpolación por polinomios de Lagrange Sea una función f (x) que conocemos en un conjunto ﬁnito de valores {xi }i=0,..,N . Es decir, sabemos que f (xi ) = fi . El polinomio interpolador de Lagrange PN (x) de f (x) en los puntos {xi }i=0,..,N , es el único polinomio de grado menor o igual que N tal que PN (xi ) = f (xi ) ∀i = 0, .., N
PN (x) se puede expresar en término de los denominados polinomios base de Lagrange P i (x), deﬁnidos como: P i (x) = ΠN i (x − xj ) j6=
Problema 24 (2 puntos) Calcular el polinomio interpolador de Lagrange P3 (x) de la función f (x) = sen(x) en los puntos 0, π , π y 3π . 2 2
estos polinomios base tienen la propiedad fundamental siguiente ½ 1 si i = j P i (xj ) = 0 si i 6= j Por tanto, el polinomio interpolador de Lagrange puede expresarse como PN (x) =
Teorema 9 El polinomio interpolador de Lagrange es el único polinomio de grado igual o inferior a N tal que PN (xi ) = f (xi ) ∀i = 0, .., N
Demostración Sea P (x) un polinomio de grado inferior o igual a N que veriﬁque que P (xi ) = f (xi ) ∀i = 0, .., N. Entonces, el polinomio Q(x) = P (x) − PN (x) es un polinomio de grado inferior o igual a N que veriﬁca que Q(xi ) = 0 y, por tanto, posee N + 1 raíces, lo cual es imposible, salvo que Q(x) sea identicamente igual a cero. Por tanto Q(x) ≡ 0 y P (x) = PN (x).
Error de interpolación de Lagrange y polinomios de Chebychev Evidentemente, al aproximar f (x) por el polinomio interpolador PN (x) en un intervalo [a, b] se comete, en general, un error de interpolación, que viene determinado por el siguiente teorema. Teorema 10 Sea f (x) una función, y PN (x) su polinomio interpolador de Lagrange en los puntos {xi }i=0,..,N ⊂ [a, b] y x ∈ [a, b], entonces f (x) − PN (x) = f N +1) (ξ) N Π (x − xi ) (N + 1)! i=0
Demostración La demostración para el intervalo [−1, 1] se encuentra en [Ki-Ch] Pg. 292-294. La demostración para un intervalo cualquiera [a, b] se obtiene fácilmente transformando el intervalo [−1, 1] en [a, b]. Por tanto, utilizando este resultado, el error de interpolación máximo viene determinado por: | f (x) − PN (x) |≤ maxx∈[a,b] f N +1) (ξ) (N + 1)!2N µ b−a 2 ¶N+1
para cualquier otra elección posible de valores de interpolación xj . e
donde ξ es un valor intermedio perteneciente a [a, b]. Demostración Si x = xi , el error de interpolación es cero y por tanto la fórmula anterior es válida. Consideremos ahora x distinto a los xi y deﬁnamos w(t) = ΠN (t − xi ) i=0 f (x) − PN (x) λ = w(x) φ(t) = f (t) − PN (t) − λw(t) La función φ(t) tiene al menos n + 1 ceros en los puntos xi y en el punto x. Por tanto, su función derivada φ0 (t) tiene al menos n ceros repartidos entre los ceros de φ(t). Análogamente, φ00 (t) tiene al menos n − 1 ceros y así sucesivamente hasta llegar a φN+1 (t), que tiene al menos 1 cero. Si llamamos ξ a dicho cero, obtenemos φN+1 (ξ) = f N+1) (ξ) − λ(N + 1)! de donde, despejando y sustituyendo λ por su valor, obtenemos el resultado del Teorema. Problema 25 (2 puntos) Calcular la expresión del error de interpolación al aproximar la función f (x) = sen(x) en el intervalo [0, 2π] interpolando en los puntos 0, π , π y 2 3π 2 , y acotarlo superiormente. La cuestión que vamos a abordar en este apartado es, en el caso en que queramos interpolar una función en un intervalo [a, b], y que nosotros podamos elegir los valores de interpolación xi , cómo elegirlos de tal forma que el error de interpolación sea mínimo. Para ello, elegiremos los puntos xi tales que ΠN (x − xi ) sea lo más pequeño posible en i=0 [a, b]. Teorema 11 Sea N ≥ 0, y un intervalo [a, b] Se consideran los puntos xi dados por µ µ ¶¶ b−a 2i + 1 1 + cos π i = 0, ..., N xi = a + 2 2N + 2
Ejemplo 6 Se considera [a, b] = [0, 1] y N = 5 (es decir 6 puntos de interpolación). Los puntos de interpolación dados por el teorema anterior son: x0 x1 x2 x3 x4 x5 = = = = = = . 982 96 . 853 55 . 629 41 . 370 59 . 146 45 1. 703 7 × 10−2
Problema 26 (2 puntos) Calcular el error máximo de interpolación en el intervalo [0, 1] al interpolar la función cos(x) en los puntos descritos en el ejemplo anterior. En el caso de que [a, b] = [−1, 1], los valores óptimos de interpolación xi dados por la fórmula anterior son las raíces de los denominados polinomios de Chebychev, TN (x), construidos de la manera siguiente: T0 (x) = 1 T1 (x) = x TN (x) = 2xTN−1 (x) − TN −2 (x) Método de diferencias de Newton para el cálculo del polinomio interpolador de Lagrange Numéricamente, el cálculo de PN (x) a través de los polinomios base necesita de la evaluación de N + 1 polinomios de grado N. Además, si queremos añadir un nuevo punto de interpolación, debemos cambiar todos los polinomios base de Lagrange. Un método más directo para el cálculo de PN (x) es el denominado método de diferencias de
Newton. El método consiste en ir calculando progresivamente los polinomios Pk (x) que interpolan la función en los puntos x0 , ..., xk de la siguiente forma: P0 (x) = a0 P1 (x) = P0 (x) + a1 (x − x0 ) P2 (x) = P1 (x) + a2 (x − x0 )(x − x1 ) ... PN (x) = PN −1 (x) + aN (x − x0 )(x − x1 )...(x − xN−1 ) A los coeﬁcientes ak los denotamos por ak = f [x0 , ..., xk ] Ejemplo 7 Vamos a interpolar la función f (x) = ex en los puntos x0 = 0, x1 = 1 y x2 = 2. P0 (x) = 1 P1 (x) = 1 + a1 x Como P1 (1) debe ser igual a e, despejando obtenemos a1 = e − 1 Por último P2 (x) = P1 (x) + a2 x(x − 1)
Demostración En primer lugar, observamos que f [xi , ..., xi+k ] indica, para cada Pk (x), el coeﬁciente que acompaña a la potencia xk en el polinomio interpolador Pk (x) para los puntos xi , ..., xi+k . Como el polinomio interpolador es único, f [xi , ..., xi+k ] no depende del orden en que tomemos los puntos xi , ..., xi+k y, por tanto: f [xi , ....., xi+k ] = f [xi+k , ....., xi ] Consideremos ahora el polinomio interpolador Qk (x) que interpola en los puntos xi+k , ..., xi , es decir, cambiando el orden de los puntos. Qk (x) se puede escribir como Qk (x) = b0 + b1 (x − xi+k ) + b2 (x − xk+i )(x − xk+i−1 ) + ... donde bj = f [xi+k , .., xi+k−j ] Por la unicidad del polinomio interpolador obtenemos que Pk (x) = Qk (x) y, por tanto ak = f [xi , ....., xi+k ] = f [xi+k , ....., xi ] = bk De nuevo, por la unicidad del polinomio interpolador, los coeﬁcientes que acompañan a la potencia xk−1 en ambos polinomios coinciden y, por tanto: ak−1 − ak
Como P2 (2) debe ser igual a e2 , despejando obtenemos e2 − P1 (2) a2 = 2 Por tanto, el polinomio P2 (x) lo expresamos como P2 (x) = 1 + (e − 1)x + e2 − 2e + 1 x(x − 1) 2
Finalmente obtenemos el resultado del teorema, teniendo en cuenta que ak−1 bk−1 = f [xi , ...., xi+k−1 ] = f [xi+k ...., xi+1 ] = f [xi+1 ...., xi+k ]
Como veremos en el teorema siguiente, los coeﬁcientes f [x0 , ..., xk ], que se denominan diferencias divididas de Newton, veriﬁcan las siguientes propiedades: f [xi ] = f (xi ) f [xi+1 ] − f [xi ] f [xi , xi+1 ] = xi+1 − xi . f [xi+1 , .., xi+k ] − f [xi , .., xi+k−1 ] f [xi , .., xi+k ] = xi+k − xi Teorema 12 Si denotamos por ak = f [x0 , .., xk ], entonces el polinomio de interpolación de Lagrange PN (x) viene dado por PN (x) =
Ejemplo 8 Sea f (x) = ex , si interpolamos f (x) en los puntos x0 = 0, x1 = 1, x2 = 2, x3 = 3, obtenemos el polinomio interpolador de la siguiente forma: f [0, 1] = e1 − 1 f [1, 2] = e2 − e1 f [2, 3] = e3 − e2 e2 − 2e + 1 f [0, 1, 2] = 2 e3 − 2e2 + e1 f [1, 2, 3] = 2 e3 − 3e2 + 3e1 − 1 f [0, 1, 2, 3] = 6 Por tanto el polinomio interpolador de Lagrange es: P3 (x) = 1 + (e − 1) x + e2 − 2e + 1 x(x − 1) + 2 e3 − 3e2 + 3e1 − 1 x(x − 1)(x − 2) 6
donde los coeﬁcientes f [xi , ..., xk ] veriﬁcan f [xi+1 , .., xi+k ] − f [xi , .., xi+k−1 ] f [xi , .., xi+k ] = xi+k − xi
En la siguiente gráﬁca se muestra la diferencia ex − P3 (x) en el intervalo [0, 3] :
0.1 0.05 0 0 -0.05 -0.1 -0.15 -0.2 -0.25 0.5 1 1.5 2 2.5 x 3
3 Problema 27 (2 puntos) Interpolar la función f (x) = 10 x2 +1 en los puntos x0 = −2, x1 = −1, x2 = 1, x3 = 2 utilizando las diferencias de Newton y evaluar el polinomio en x = 0 utilizando el algoritmo de Horner. Problema 28 (2 puntos) Calcular el polinomio interpolador de Lagrange P3 (x) de la función f (x) = sen(x) en los puntos 0, π , π y 3π utilizando las diferencias divididas 2 2 de Newton.
PARAMETER(Nmax=1000) DIMENSION A(0:1000),F(0:1000),X(0:1000) PRINT *,’Introducir No. Ptos Interp.’ READ *,N N=N-1 PRINT *,’Introducir Ptos Interpol.’ DO 1 K=0,N READ *, X(K) PRINT *,’Introducir Valores de F()’ DO 2 K=0,N READ *,F(K) IF(IDIFNEWTON(A,X,F,N).EQ.1) THEN PRINT *,’Puntos de Interpolacion repetidos’ STOP ENDIF PRINT *,’Coef. Polinomio’ DO 3 K=0,N PRINT *,A(K) PRINT *,’Test de Comprobacion’ DO 4 K=0,N PRINT *,X(K),F(K),EVDIFNEWTON(A,X,X(K),N) END FUNCTION IDIFNEWTON(A,X,F,N) Parameter(Nmax=1000) DIMENSION A(0:*),X(0:*),F(0:*),B(0:Nmax) DO 1 K=0,N B(K)=F(K) A(0)=F(0) DO 2 K=1,N DO 3 L=0,N-K IF (X(K+L).EQ.X(L)) THEN IDIFNEWTON=1 RETURN ENDIF B(L)=(B(L+1)-B(L))/(X(K+L)-X(L)) CONTINUE A(K)=B(0) CONTINUE IDIFNEWTON=0 END
Problema 29 (3 puntos) Calcular el polinomio interpolador de Lagrange P3 (x) de la función f (x) = 2x en los puntos 0, 1, 3 y 4 utilizando las diferencias divididas de Newton. Expresar el polinomio tomando en primer lugar x0 = 0, x1 = 1, x2 = 3 y x3 = 4 y, en segundo lugar, x0 = 4, x1 = 3, x2 = 1 y x3 = 0. Problema 30 (3 puntos) Dada una función f (x) y una secuencia de valores xn , aproximar f (x) por la parábola que pasa por los puntos (xn−1 , f (xn−1 )) , (xn−2 , f (xn−2 )) y (xn−3 , f (xn−3 )). Calcular posteriormente las derivadas del polinomio y comprobar que coinciden con las fórmulas dadas en el método de Müller para el cálculo de las derivadas f 00 (xn−1 ) y f 0 (xn−1 ). Programa 10 Programa en Fortran 77 donde se deﬁnen las funciones IDIFNEWTON, que a partir del vector X(0 : N ) de puntos de interpolación y el vector F (0 : N ) de valores de la función f (x) en los puntos de interpolación, devuelve el vector A(0 : N ) de coeﬁcientes de diferencias divididas que deﬁnen el polinomio de Lagrange (A(K) = f [x0 , x1 , .., xK ] ), y la función EVDIFNEWTON(A,X,X0,N) que a partir de los coeﬁcientes dados por el vector A(0 : N ) y el conjunto de puntos de interpolación, devuelve el valor de la evaluación del polinomio de Lagrange en el punto X0.
FUNCTION EVDIFNEWTON(A,X,X0,N) DIMENSION A(0:*),X(0:*) EVDIFNEWTON=A(N) DO 1 K=N-1,0,-1 1 EVDIFNEWTON=EVDIFNEWTON*(X0X(K))+A(K) END
Implementación de funciones elementales Una vez deﬁnida una aritmética en precisión ﬁnita y las 4 operaciones básicas (suma, resta, multiplicación, división), es necesario deﬁnir, a partir de estas operaciones, las funciones elementales que todos usamos, como son: la raíz √ cuadrada x, las funciones trigonométricas: sen(x) cos(x) y tan(x), la función ln(x), la función ex , la función xy , etc. Las técnicas elementales para deﬁnir estas funciones consisten en utilizar la interpolación polinómica, los desarrollos de Taylor y los algoritmos de búsqueda de ceros (como √ vimos anteriormente para x). Aproximación de la exponencial ex Un número real x siempre se puede expresar como x = m + x0 , donde m es un número entero y x0 ∈ [0, 1]. Dado que 0 ex = em ex podemos descomponer el cálculo de ex en el cálculo, por un lado, de em , donde al ser m un entero el cálculo es inmediato a partir de multiplicaciones sucesivas de potencias naturales de e ó e−1 (si m<0), y por otro, en el cálculo 0 de ex para x0 ∈ [0, 1]. Utilizando como puntos de interpolación los asociados a los polinomios de Chebychev: µ µ ¶¶ 1 2i + 1 xi = 1 + cos π i = 0, ..., N 2 2N + 2 obtenemos que el error relativo veriﬁca que: µ ¶N+1 0 e 1 | ex − PN (x) | ≤ ex0 (N + 1)!2N 2 Para N = 6, el error relativo es menor que 6.6 × 10−8 y, por tanto, del mismo orden que la unidad de redondeo u en una aritmética de 32 bits. Asi, tomando un polinomio de grado N = 6, es decir 7 puntos de interpolación, obtenemos ya la mejor aproximación posible de ex en el intervalo [0, 1] en una aritmética de 32 bits. Práctica 3 (Aproximación de ex , 2 horas) Crear una función en Fortran 77 que devuelva el valor de ex con x ∈ [0, 1] utilizando el polinomio de Lagrange P6 (x) que interpola a ex en los puntos: 1 xi = 2 µ µ ¶¶ 2i + 1 1 + cos π 14
Aproximación de funciones trigonométricas Utilizaremos como modelo las funciones f (x) = cos(x) y f (x) = sen(x). Puesto que estas funciones son 2π periódicas, utilizando algunas relaciones trigonométricas es suﬁciente deﬁnir las funciones cos(x) y sen(x) en el intervalo [0, π ] y a partir de ellas deﬁnir las funciones para 4 cualquier valor x (en radianes). Efectivamente, denotemos por cos[0, π ] (x) y sen[0, π ] (x) a las funciones trigonométri4 4 cas deﬁnidas sobre el intervalo [0, π ]. Podemos deﬁnir en4 tonces las siguientes funciones: ½ cos[0, π ] (x) si x ≤ π 4 4 cos[0, π ] (x) = π 2 sen[0, π ] ( 2 − x) si x > π 4 4 sen[0, π ] (x) = 2 cos[0,π] (x) = ½ sen[0, π ] (x) si x ≤ 4 π π ( cos[0, 4 ] 2 − x) si x >
sen[0,π] (x) = cos[0,2π] (x) = sen[0,2π] (x) =
cos[0, π ] (x) si x ≤ 2 − cos[0, π ] (π − x) si x > 2 sen[0, π ] (x) si x ≤ 2 sen[0, π ] (π − x) si x > 2
cos[0,π] (x) si x ≤ π cos[0,π] (2π − x) si x > π sen[0,π] (x) si x ≤ π −sen[0,π] (2π − x) si x > π
El desarrollo en Serie de Taylor centrado en 0 del cos(x) es: cos(x) u Pn (x) = 1. − x2 x4 x2n + + .... + (−1)n 2 4! (2n)!
y el error máximo cometido por el desarrollo de Taylor en un punto x ∈ [0, π ] es 4 | Pn (x) − cos(x) |≤ sen(x) (x)2n+1 (2n + 1)!
La ventaja de utilizar el desarrollo de Taylor centrado en 0 es que las potencias impares de x no aparecen, lo que simpliﬁca el cálculo numérico. El error relativo es | Pn (x) − cos(x) | (x)2n+1 ≤ tan(x) cos(x) (2n + 1)! Además, como tan(x) es creciente en [0, π ], el valor 4 máximo del error se encuentra en x = π . Por ejemplo, para 4 n = 5 obtenemos que el error relativo máximo cometido en x = π es del orden de 4 ¡ π ¢2∗5+1 π tan( ) 4 = 1. 8 × 10−9 4 (2 ∗ 5 + 1)! Por tanto, si trabajamos con una aritmética de 32 bits, cuya unidad de redondeo u es del orden de 10−8 , tenemos que con n = 5 obtenemos una aproximación del cos(x) que es la mejor posible dentro de esta aritmética y no tendría sentido aumentar el valor de n.
i = 0, ..., 6
Comprobar que el polinomio esta bien construido, es decir que P6 (xi ) = exi para todos los xi . Introducir por teclado un valor x, evaluar y mostrar P6 (x) y ex . Utilizar x = 0, 1, 0.5, 2 y 3. Nota: Utilizar las funciones de an.h IDIFNEWTON(.), que calcula el polinomio interpolador a partir de los puntos y valores de interpolación, y EVDIFNEWTON(.), que evalua el polinomio interpolador en un punto.
Problema 31 (3 puntos) Aproximar la función sen(x) en el intervalo [0, π ] utilizando el desarrollo de Taylor y 4 calcular el valor de n a partir del cual la aproximación es la mejor posible dentro de una aritmética de 32 bits.
Para N = 10 el error máximo es 3. 973 6 × 10−8 , que es menor que la unidad de redondeo u y, por tanto, en una aritmética de 32 bits tendríamos la mejor aproximación posible de la función ln(x). Problema 34 (1 punto) ¿Cómo se puede obtener la función y x , donde x e y son números reales, utilizando las funciones ex y ln(x)?
Problema 32 (2 puntos) Demostrar que, utilizando relaciones trigonométricas, es posible calcular las funciones sen(x) y cos(x) para cualquier x (en radianes), utilizando únicamente su valor en el intervalo [0, π ]. 8
Problema 33 (3 puntos) Calcular los polinomios necesarios para interpolar las funciones trigonométricas cos(x) y sen(x) en el intervalo [0, π ] en una aritmética de 32 bits. 8
ANÁLISIS NUMÉRICO MATRICIAL I En esta primera sección dedicada a la resolución de sistemas de ecuaciones lineales, estudiaremos los métodos directos clásicos para la resolución de un sistema de ecuaciones de la forma Au = b donde A = (ai,j ) es una matriz de N xN, b = (bi ) es un vector de tamaño N que determina los términos independientes, y u = (ui ) es el vector solución buscado. Método de Gauss Este método, aunque no es de los más rápidos, tiene la gran ventaja de que se puede aplicar a todo tipo de matrices, algo que, como veremos en el futuro, no ocurre con otros métodos más rápidos, pero que requieren, por ejemplo, que la matriz sea simétrica o deﬁnida positiva. El método de Gauss se basa en transformar el sistema Au = b en un sistema equivalente A0 u = b0 tal que la solución sea la misma y que la matriz A0 sea triangular superior, es decir, que tenga valores nulos de la diagonal hacia abajo. Una vez obtenidos la matriz A0 y el vector b0 , el cálculo de la solución u es inmediata, siguiendo un remonte de las variables a través del siguiente esquema recursivo: uN = bk −
Aproximación de la función ln(x) Como hemos visto anteriormente, un número x real en una aritmética de precisión ﬁnita viene expresado habitualmente como Ã t ! X an m x=2 2n n=1
n donde m es un número entero, a1 = 1 y para´ > 1 an = 0 ³P t an ó an = 1. Por tanto, el número es mayor o n=1 2n igual que 1 y menor que 1. Aplicando las propiedades del 2 ln(x) obtenemos que Ã t ! X an ln(x) = m ln(2) + ln 2n n=1
Dado que el número ln(2) es una constante que supondremos calculada anteriormente ( ln(2) ∼ .6931471806), = podemos reducir el cálculo del ln(x) al rango de valores 1 2 ≤ x ≤ 1. Utilizaremos los puntos de interpolación generados por los polinomios de Chebychev, que para el intervalo [ 1 , 1] son: 2 µ µ ¶¶ 1 1 2i + 1 xi = + 1 + cos π i = 0, ..., N 2 4 2N + 2 Dado que ln(1) = 0, para minimizar el error relativo añadiremos como punto interpolante xN +1 = 1. El error de interpolación relativo entre PN +1 (x) y ln(x) es: | (x − 1)ΠN (x − xi ) | | ln(x) − PN+1 (x) | i=0 = N +1 | ln(x) | ξ (N + 2) | ln(x) | donde ξ ∈ [ 1 , 1]. Además se tiene que en el intervalo [ 1 , 1] 2 2 |x−1| ≤1 | ln(x) | Por tanto: | ln(x) − PN +1 (x) | 2N+1 ≤ | ln(x) | (N + 2) µ ¶N +1 1 1 4 2N
b0 N a0 N,N k = N − 1, .., 1
l=k+1 a0 k,k
a0 ul k,l
Problema 35 (2 puntos) Calcular el número de operaciones básicas (sumas, restas, multiplicaciones y divisiones) necesarias para realizar un remonte como el presentado arriba en función de la dimensión N . Para obtener A0 y b0 se calcula, en primer lugar, el valor máximo en valor absoluto de la primera columna de A, denominado pivote. A continuación, se intercambia la primera ﬁla de A con la ﬁla donde se encuentra el pivote, y se hace lo mismo con el vector b, para que el sistema sea equivalente. A continuación, se multiplica la primera ﬁla de A por el valor −ak y se suma a la ﬁla k − esima ´ a1 de A para k = 2, ..., N. Se hace lo mismo para el vector
b, y con ello habremos obtenido un sistema equivalente tal que la primera columna es cero de la diagonal hacia abajo. Volvemos ahora a hacer lo mismo para convertir la segunda columna cero de la diagonal para abajo, y así sucesivamente hasta llegar a la mencionada matriz A0 .
Ejemplo 9 Ejemplo de descomposición según el método de Gauss. Se considera el sistema ⎛ ⎞⎛ ⎞ ⎛ ⎞ −2 −2 0 u1 0 ⎝ 6 18 12 ⎠ ⎝ u2 ⎠ = ⎝ 24 ⎠ 3 11 7 8 u3
de la misma forma, forma siguiente ⎛ ⎞ ⎛ 0 ⎝ 24 ⎠ → ⎝ 8
Programa 11 Programa en fortran 77 que implementa el método de Gauss. Se deﬁne una función IGAU SS que tiene como parámetros la matriz A, el vector independiente b, un vector auxiliar N row, la dimensión del sistema y la dimensión máxima admitida. La función devuelve un valor entero M que indica si se ha terminado correctamente (M = 0) o incorrectamente (M = 1, 2). En el caso en que se ha terminado correctamente, la solución se devuelve en el propio vector b. Parameter(Nmax=1000) DIMENSION A(Nmax,Nmax),B(Nmax),Nrow(Nmax) PRINT *, ’Introducir Dimension’ READ *,N PRINT *,’Introducir matriz’ DO 2 K=1,N DO 2 L=1,N READ *,A(K,L) PRINT *,’Introducir vector’ DO 3 K=1,N READ *,B(K) M=IGAUSS(A,B,N,Nrow,Nmax) IF(M.EQ.0) THEN PRINT *,’SOLUCION’ DO 1 K=1,N PRINT *,B(K) CONTINUE ELSE IF(M.EQ.1)THEN PRINT *,’Una columna es toda cero’ ELSE PRINT *,’A(N,N)=0’ ENDIF END
y el remonte da como solución u1 = 6, u2 = −6, u3 = 8. 1
Problema 36 (2 puntos) Resolver por el método de Gauss el sistema¶ µ µ ¶ µ ¶ −1 2 x 3 = 2 −1 y 0 Problema 37 (3 puntos) Calcular el número de operaciones básicas necesarias para descomponer el sistema Au = b en el sistema A0 u = b0 utilizando el método de Gauss, teniendo en cuenta la siguiente relación:
Problema 38 (2 puntos) Implementar en FORTRAN la funcion IDESCEN SO(A, b, u, N, N max) que resuelve un sistema, donde A es una matriz triangular inferior, b es el vector de términos independientes, u el vector solución, N es la dimensión real del sistema y N max la dimensión que se utilizó para reservar la memoria de la matriz A. La función devuelve 0 si termina correctamente y 1 en caso contrario. Nota Importante: Las líneas de código tienen que ir todas numeradas y no pueden superar las 15 líneas.
FUNCTION IGAUSS(A,B,N,Nrow,Nmax) PARAMETER (NmaxGAUSS=1000) DIMENSION A(Nmax,*),B(*),Nrow(*) DIMENSION U(NmaxGAUSS) IF (N.GT .NmaxGAUSS) THEN IGAUSS=3 PRINT *,’DIMENSION DEL SISTEMA MAYOR DE LA PERMITIDA’ RETURN ENDIF DO 1 K=1,N 1 Nrow(K)=K DO 2 K=1,N-1 XMax=ABS(A(Nrow(K),K)) M=K DO 3 L=K+1,N
IF(ABS(A(Nrow(L),K)).GT.XMax) THEN XMax=ABS(A(Nrow(L),K)) M=L ENDIF 3 CONTINUE IF(XMax.LT.(2.**(-100))) THEN IGAUSS=1 RETURN ENDIF IF(K.NE.M) THEN MP=Nrow(K) Nrow(K)=Nrow(M) Nrow(M)=MP ENDIF DO 4 L=K+1,N C=A(Nrow(L),K)/A(Nrow(K),K) DO 5 M=K,N A(Nrow(L),M)=A(Nrow(L),M)C*A(Nrow(K),M) 5 CONTINUE B(Nrow(L))=B(Nrow(L))-C*B(Nrow(K)) 4 CONTINUE 2 CONTINUE IF(ABS(A(Nrow(N),N)).LT.(2.**(-100))) THEN IGAUSS=2 RETURN ENDIF U(N)=B(Nrow(N))/A(Nrow(N),N) DO 6 K=N-1,1,-1 C=0 DO 7 L=K+1,N C=C+A(Nrow(K),L)*U(L) 7 CONTINUE U(K)=(B(Nrow(K))-C)/A(Nrow(K),K) 6 CONTINUE DO I=1,N B(I)=U(I) ENDDO IGAUSS=0 END
donde N es la dimensión del sistema y ErrorSistema representa el error relativo medio al resolver el sistema. En el denominador se añade 1 para evitar las posibles divisiones por 0. Cuanto más pequeño sea ErrorSistema, mejor aproximada estará la solución del sistema.
Método de Cholesky Este método sólo se puede aplicar a matrices simétricas y deﬁnidas positivas. El siguiente teorema da 3 posibles deﬁniciones equivalentes de una matriz deﬁnida positiva. Teorema 13 Sea A una matriz simétrica, las 3 siguientes aﬁrmaciones son deﬁniciones equivalentes a que una matriz sea deﬁnida positiva (i) ∀ v ∈ <N − {0} se cumple que t vAv > 0. (ii) Todos los autovalores λ de A son positivos. (iii) Los determinantes de todos los menores principales de A son positivos. El método de Cholesky se basa en descomponer la matriz A en la forma: A = B·B t donde B es una matriz ⎛ b1,1 ⎜ b2,1 ⎜ B = ⎜ b3,1 ⎜ ⎝ . bn,1 triangular inferior. 0 b2,2 b3,2 . bn,2 0 0 b3,3 . bn,3 . 0 . . . . . 0 . bn,n ⎞ ⎟ ⎟ ⎟ ⎟ ⎠
Problema 40 (2 puntos) Demostrar que si A = B · B t (B triangular inferior) y |B| 6= 0, entonces A es simétrica y deﬁnida positiva. Problema 41 (2 puntos) Descomponer la siguiente matriz A por el método de Cholesky: ⎛ ⎞ 1 1 4 A=⎝ 1 5 6 ⎠ 4 6 26 De forma general, el algoritmo para calcular B es el siguiente Para i =r ..., N 1, ³ ´ P bi,i = ai,i − i−1 b2 k=1 i,k
Estimación del error de un método para resolver sistemas Para estimar la ﬁabilidad de la solución numérica de un sistema de ecuaciones, haremos lo siguiente: dada una matriz A, un vector de términos independientes b y un vector solución u, calculado utilizando alguna técnica numérica, si la solución es perfecta entonces Au − b = 0. Ahora bien, esto no suele suceder, porque los errores de redondeo y de cálculo producen que esta estimación no sea exacta. Para estimar el error cometido al resolver el sistema utilizaremos la expresión siguiente, donde e es el vector e = Au − b : ErrorSistema = 1 X |ei | N |bi | + 1
Para j = i + 1, ..., N ³ ´ P 1 bj,i = bi,i aj,i − i−1 bj,k bi,k k=1 Fin Para j Fin Para i
El interés de descomponer una matriz A por el método de Cholesky es que, a continuación, es muy sencillo resolver el sistema de ecuaciones Au = b. Efectivamente, basta descomponer el sistema de la siguiente forma: Bz = b Bu = z
Ambos sistemas se resuelvan rápidamente haciendo un remonte y un descenso. Nota: Normalmente, para evitar tener que almacenar dos matrices, una para B y otra para B t , se almacena todo en una única matriz B, simétrica, escribiendo en la parte triangular superior de B la parte correspondiente a B t . Problema 42 (2 puntos) Calcular el número de operaciones necesarias para resolver un sistema por el método de Cholesky.
Problema 43 (2 puntos) Demostrar que a partir de un método para resolver sistemas de ecuaciones se puede construir de forma inmediata un método para calcular la inversa A−1 de una matriz A.
Hay que hacer una versión en simple precisión y otra versión en doble precisión donde todas las variables que empiecen por D sean de doble precisión (la matriz A y el vector B siempre serán de simple precisión). Hay que resolver los sistemas ejemplo y calcular el ErrorSistema, tanto en doble como en simple precisión. Nota: El programa debe permitir introducir el sistema directamente por teclado o desde disco duro, utilizando las funciones deﬁnidas en an.h. Resolver los siguientes sistemas ejemplo: ⎞⎛ ⎞ ⎛ ⎞ ⎛ x 6 1 1 4 1. ⎝ 1 5 6 ⎠ ⎝ y ⎠ = ⎝ 12 ⎠ 4 6 26 z 36 ⎛ ⎞⎛ ⎞ ⎛ ⎞ 1 1 4 x 6 2. ⎝ 1 1 4 ⎠ ⎝ y ⎠ = ⎝ 12 ⎠ 4 4 17 z 36 ⎛ ⎞⎛ ⎞ ⎛ ⎞ 1 1 4 x 6 3. ⎝ 1 5 6 ⎠ ⎝ y ⎠ = ⎝ 12 ⎠ 4 6 17 z 36 ⎛ ⎞⎛ ⎞ ⎛ ⎞ 2 −1 0 x 1 4. ⎝ −1 2 −1 ⎠ ⎝ y ⎠ = ⎝ 0 ⎠ 0 −1 −4 z −5 Resolver también los sistemas que aparecen en el directorio /users/asignaturas/ii-an de la máquina serdis.dis.ulpgc.es. En este directorio hay tres ejemplos de sistemas de dimensión 10, 100 y 500. Los 3 sistemas corresponden a matrices simétricas y deﬁnidas positivas. Estos archivos ejemplo sólo se pueden utilizar al compilar el programa en serdis bajo UNIX. Si utilizamos Linux, el programa no reconocerá el formato de los archivos.
Práctica 4 (Método de Cholesky, 6 horas) Implementar en Fortran 77 las siguientes funciones : • FUNCTION ICHOLESKY_FACTORIZACION (A,DB,N,Nmax): Calcula la descomposición de Cholesky de A y la devuelve en la matriz DB. Devuelve 0 si termina bien y −1 en caso contrario. • FUNCTION IDESCENSO(DB,DZ,B,N,Nmax): Resuelve un sistema triangular inferior, donde DB es la matriz, B es el término independiente y DZ es el vector donde devuelve la solución. Devuelve 0 si termina bien y −1 en caso contrario. • FUNCTION IREMONTE(DB,DU,DZ,N,Nmax): Resuelve un sistema triangular superior, donde DB es la matriz, DZ es el término independiente y DU es el vector donde devuelve la solución. Devuelve 0 si termina bien y −1 en caso contrario • FUNCTION ERROR_SISTEMA(A,DU,B,N,Nmax): Devuelve el error cometido al resolver el sistema dado por la expresión ErrorSistema de la sección anterior. • FUNCTION ICHOLESKY(A,B,DU,N,Nmax): Resuelve un sistema por el método de Cholesky y devuelve la solución en DU. Devuelve 0 si termina bien y −1 en caso contrario.
Método de Crout para matrices tridiagonales El caso de sistemas de ecuaciones con matrices A tridiagonales posee una forma especialmente simple de factorización. Vamos a descomponer A en el producto de dos matrices triangulares de la forma siguiente: ⎛ ⎞ a1 b1 . 0 ⎜ c1 a2 . 0 ⎟ ⎜ ⎟= ⎝ 0 . . bN −1 ⎠ 0 . cN −1 aN ⎞ ⎞⎛ ⎛ 1 u1 . . 0 0 l1 0 ⎜ m1 l2 0 ⎟ . 0 ⎟⎜ 0 1 . ⎟ ⎟⎜ ⎜ ⎝ 0 . . 0 ⎠ ⎝ 0 . . uN −1 ⎠ 0 . mN−1 lN 0 . 0 1
Los vectores mi , li, y ui se calculan utilizando el esquema: l1 = a1 1 u1 = b1 l Para i = 2, .., N − 1
mi−1 = ci−1 li = ai − mi−1 ui−1 i ui = bi l Fin Para mN−1 = cN −1 lN = aN − mN −1 uN −1 Problema 44 (3 puntos) Demostrar el algoritmo de Crout para descomponer matrices tridiagonales.
DIMENSION Vector(*) OPEN(1,FILE=String, form=’UNFORMATTED’,ERR=1) GOTO 2 1 OPEN(1,FILE=String, form=’UNFORMATTED’) 2 A=Ndimension WRITE(1) A DO 3 K=1,Ndimension WRITE(1) Vector(K) 3 CONTINUE CLOSE(1) END
STATUS=’NEW’,
STATUS=’OLD’,
Subrutinas en Fortran 77 para la lectura y escritura en disco de vectores y matrices Programa 12 Subrutinas en Fortran 77 que leen y escriben en disco vectores y matrices. ****** Función que lee un vector de disco de nombre String, lo almacena ****** en la tabla Vector y devuelve la dimension del vector FUNCTION LeerVector(String,Vector) CHARACTER * (*) String DIMENSION Vector(*) OPEN(1,FILE=String,STATUS=’OLD’, form=’UNFORMATTED’) READ(1) A Ndimension=INT(A) PRINT *,’DIMENSION DEL VECTOR ’,String,’=’,Ndimension DO 1 K=1,Ndimension READ(1) Vector(K) 1 CONTINUE CLOSE(1) LeerVector=Ndimension END
****** Función que lee una matriz cuadrada de disco de nombre String, lo almacena ****** en la tabla A y devuelve la dimension de la matriz FUNCTION LeerMatriz(String,A,Nmax) CHARACTER * (*) String DIMENSION A(Nmax,*) OPEN(1,FILE=String, STATUS=’OLD’, form=’UNFORMATTED’) READ(1) B Ndimension=INT(B) PRINT *,’DIMENSION DE LA MATRIZ ’, String,’=’, Ndimension DO 1 K=1,Ndimension DO 2 L=1,Ndimension READ(1) A(K,L) 2 CONTINUE 1 CONTINUE CLOSE(1) LeerMatriz=Ndimension END
***** Procedimiento que escribe en el ﬁchero String el vector Vector ***** de dimension Ndimension SUBROUTINE EscribirVector(String,Vector,Ndimension) CHARACTER * (*) String
***** Procedimiento que escribe en el ﬁchero String la matriz A ***** cuadrada de dimension Ndimension SUBROUTINE EscribirMatriz(String,A,Ndimension,Nmax) CHARACTER * (*) String DIMENSION A(Nmax,*) OPEN(1,FILE=String, STATUS=’NEW’, form=’UNFORMATTED’, ERR=1) GOTO 2 1 OPEN(1, FILE=String, STATUS=’OLD’, form=’UNFORMATTED’) 2 B=Ndimension WRITE(1) B DO 3 K=1,Ndimension DO 4 L=1,Ndimension WRITE(1) A(K,L) 4 CONTINUE 3 CONTINUE CLOSE(1)
DIFERENCIACIÓN E INTEGRACIÓN NUMÉRICA Una fórmula de diferenciación numérica es un procedimiento que permite aproximar la derivada de la función f (x) en un punto xi. utilizando el valor de f (x) en otros puntos vecinos a xi . Por otro lado, una fórmula de integración numérica es un procedimiento que permite aproximar el valor de la integral en un intervalo [a, b] a partir de la evaluación de f (x) en algunos puntos incluidos en el intervalo [a, b]. Diferenciación Numérica La manera habitual de aproximar la derivada de una función f (x) en un punto xi consiste en utilizar el desarrollo de Taylor centrado en xi : f (x) = f (xi ) + f 0 (xi ) f N) (xi ) (x − xi ) + ... + (x − xi )N + ... 1! N!
Programa 13 Programa en Fortran 77 donde se describe un ejemplo de lectura/escritura de vectores y matrices. INCLUDE ’an.h’ PARAMETER(Nmax=100) CHARACTER * 10,String DIMENSION V(Nmax), A(Nmax,Nmax) String=’vector.dat’ Ndimension=LeerVector(String,V) PRINT *,Ndimension DO 1 K=1,Ndimension PRINT *,V(K) V(K)=2*V(K) CONTINUE CALL EscribirVector(’vector3.dat’,V,Ndimension) N=LeerVector(’vector3.dat’,V) PRINT *,’N=’,N DO 2 k=1,N PRINT *,V(K) CONTINUE Ndimension=3 DO 3 K=1,Ndimension DO 4 L=1,Ndimension A(L,K)=L+K PRINT *,A(L,K) CONTINUE PRINT * CONTINUE CALL EscribirMatriz(’matriz.dat’,A,Ndimension,Nmax) Ndimension=LeerMatriz(’matriz.dat’,A,Nmax) PRINT *,Ndimension DO 5 K=1,Ndimension DO 6 L=1,Ndimension PRINT *,A(L,K) CONTINUE PRINT * CONTINUE END
Si tomamos un punto xj 6= xi , truncamos el desarrollo de Taylor y despejamos, obtenemos la siguiente expresión: f 0 (xi ) ≈ f (xj ) − f (xi ) + O (|xj − xi |) xj − xi
donde O (|xj − xi |) indica, básicamente, que el error cometido es una suma de potencias de |xj − xi | en la que la potencia más pequeña es 1. Se denomina orden de la aproximación a la potencia más pequeña que aparece en el término del error. Por lo tanto, en este caso, diremos que el orden de aproximación es 1. Si xj > xi , entonces la derivada se calcula hacia adelante, mientras que si xj < xi , la derivada se calcula hacia atrás. Ejemplo 10 Veremos en este ejemplo como, cuanto más próximo esté el punto xj al punto xi , mejor será el valor aproximado de la derivada. Consideremos la función f (x) = x3 . La derivada de f (x) en x = 1 es f 0 (1) = 3 Si tomamos xi = 1 y xj = 2 en la fórmula anterior, obtenemos la aproximación f 0 (1) ≈ 23 − 13 =7 2−1
Si tomamos ahora xj = 1.1, obtenemos f 0 (1) ≈ 1.13 − 13 = 3. 31 1.1 − 1
Nota la declaración IN CLU DE 0an.h0 incluye el ﬁchero an.h que se encuentra en el directorio de trabajo, en el cuerpo del programa. La declaración CHARACT ER ∗ 10, String deﬁne un string de caracteres de tamaño 10.
que está mucho más próximo al valor real.
Problema 47 (2 puntos) Calcular analítica y numéricamente la matriz gradiente en el punto (1, 1) (utilizar h = 0.1) de la función: ½ 2 x + y2 − 1 f (x, y) = x−y
Problema 48 (3 puntos) Dados 3 puntos distintos xl , xi , xr , demostrar que la fórmula f (xi ) ≈
Problema 52 (2 puntos) Considerar en el problema anterior que xl = xi − h, y xr = xi + h. Deducir como queda la fórmula anterior para aproximar la derivada segunda, y demostrar que, en este caso, el orden de aproximación es 2. Problema 53 (3 puntos) Dados 3 puntos xl < xi < xr , calcular el polinomio de Lagrange que interpola a f (x) en esos 3 puntos, calcular la derivada segunda de ese polinomio en xi , y comprobar que da la misma fórmula que utilizando los desarrollos de Taylor. Problema 54 (2 puntos) Calcular una aproximación de la derivada primera y segunda de una función f (x) en x = 0, teniendo en cuenta que f (0) = 1, f (1) = 0, f (4) = 9
aproxima la derivada de f 0 (xi ) con un orden de aproximación de 2. Nótese que, si xr = xi + h, y xl = xi − h, entonces la fórmula anterior resulta f 0 (xi ) = f (xi + h) − f (xi − h) 2h
que es una conocida fórmula de diferencias centradas. Ejemplo 11 Veremos en este ejemplo como, utilizando la expresión anterior para aproximar la derivada de f (x) = x3 en x = 1, la precisión es mayor que con la fórmula anterior. Por ejemplo, si tomamos xi = 1 y h = 1, la expresión anterior nos da f 0 (1) ≈ 23 − 03 =4 2
Diferenciación numérica en dimensiones superiores Estudiaremos, en este apartado, la aproximación de las derivadas de una función de varias variables. Para simpliﬁcar la exposición, supondremos que la dimensión es 2. Para discretizar las derivadas de una función F (x, y), se utilizan los desarrollos de Taylor siguientes en 2 variables. (x,y) Utilizaremos la siguiente nomenclatura: Fx = ∂F∂x , Fy =
∂F (x,y) , ∂y ∂ 2 F (x,y) ∂y2
Si tomamos ahora xj = 0.1 1.13 − 0.93 f 0 (1) ≈ = 3. 01 0.2 que está más próximo al valor real que utilizando la primera fórmula. En general, comprobamos que, cuanto mayor es el orden de una fórmula de aproximación, más preciso es el valor de la derivada. Nota: Utilizar el desarrollo de Taylor para aproximar f 0 (x) es equivalente a interpolar f (x) con el polinomio de Lagrange y posteriormente derivar el polinomio. Problema 49 (3 puntos) Dados 3 puntos distintos xl , xi , xr , calcular el polinomio de Lagrange que interpola a f (x) en esos 3 puntos, calcular la derivada de ese polinomio en xi , y comprobar que da la misma fórmula que la presentada en el problema anterior. Problema 50 (2 puntos) Calcular una aproximación de la derivada tercera f 000 (xi ) de una función f (x) en un punto xi , utilizando f (xi ), f (xi + h), f (xi − h), f (xi − 2h). Problema 51 (3 puntos) Dados 3 puntos, demostrar que la fórmula f (xi ) ≈ 2
∂ 2 F (x,y) , ∂x2
∂ 2 F (x,y) ∂x∂y ,
1. F (x + h, y) = F + hFx + 2. F (x − h, y) = F − hFx + 3. F (x, y + l) = F + lFy + 4. F (x, y − l) = F − lFy +
6. F (x − h, y − l) = F − hFx − lFy + 1 (h2 Fxx + 2hlFxy + 2 ¢3 ¡ l2 Fyy ) + O( h2 + l2 2 ) 7. F (x + h, y − l) = F + hFx − lFy + 1 (h2 Fxx − 2hlFxy + 2 ¡ ¢3 l2 Fyy ) + O( h2 + l2 2 ) 8. F (x − h, y + l) = F − hFx + lFy + 1 (h2 Fxx − 2hlFxy + 2 ¡ ¢3 l2 Fyy ) + O( h2 + l2 2 ) Prestaremos particular atención a dos operadores diferenciales que se utilizan con frecuencia en la práctica: El gradiente ∇F (x, y) = (Fx (x, y), Fy (x, y)), que es el vector de derivadas parciales, y el Laplaciano ∆F (x, y) = Fxx (x, y) + Fyy (x, y). Utilizaremos la notación Fi,j ∼ = F (hi, lj).
5. F (x + h, y + l) = F + hFx + lFy + 1 (h2 Fxx + 2hlFxy + 2 ¡ ¢3 l2 Fyy ) + O( h2 + l2 2 )
aproxima la derivada segunda de f (x) en xi con un orden de aproximación de 1.
Discretización del Laplaciano Para discretizar el operador ∆F en un entorno de 3 × 3 puntos, pueden utilizarse diferentes esquemas. Para simpliﬁcar, supondremos que l = h. Problema 55 (3 puntos) Demostrar, utilizando el desarrollo de Taylor, que las siguientes expresiones son discretizaciones del laplaciano: ∆F ∆F = = Fi+1,j+1 + Fi−1,j+1 + Fi−1,j−1 + Fi+1,j−1 − 4Fi,j 2h2 Fi+1,j + Fi−1,j + Fi,j+1 + Fi,j−1 − 4Fi,j h2
Problema 56 (2 puntos) Calcular una aproximación del laplaciano de una función F (x, y) en el punto (x, y) = (0, 0) conociendo los siguientes valores: F (0, 0) = 0, F ( 1 , 0) = 1 , F (− 1 , 0) = 1 , F (0, 1 ) = 1 , F (0, − 1 ) = 1 , 2 4 2 4 2 4 2 4 F ( 1 , 1 ) = 1 , F (− 1 , − 1 ) = 1 , F (− 1 , 1 ) = 1 , F ( 1 , − 1 ) = 2 2 2 2 2 2 2 2 2 2 2 1 2. Discretización del gradiente Siguiendo el desarrollo de Taylor mostrado anteriormente, obtenemos la siguiente expresión para el gradiente: (Fi+1,j − Fi−1,j ) + 2h (Fi+1,j+1 − Fi−1,j+1 + Fi+1,j−1 − Fi−1,j−1 ) +γ 4h (Fi,j )x = (1 − γ) (Fi,j+1 − Fi,j−1 ) + 2h (Fi+1,j+1 − Fi+1,j−1 + Fi−1,j+1 − Fi−1,j−1 ) +γ 4h donde γ es, de nuevo, un parámetro a elegir. Teniendo en cuenta que la norma euclídea del gradiente es invariante por rotaciones, lo será en particular para rotaciones de 45 grados, de donde deducimos, utilizando el mismo √ argumento que para el ∆F , que γ = 2 − 2. Por lo tanto, estamos calculando Fx utilizando la máscara √ √ −(2 − 2) 0 (2 − 2) √ √ 1 −2( 2 √ 1) 0 2( 2 √ 1) − − 4h −(2 − 2) 0 (2 − 2) (Fi,j )y = (1 − γ) y Fy utilizando √ −(2 − 2) 1 0√ 4h (2 − 2) √ 2( 2 − 1) √0 −2( 2 − 1) √ −(2 − 2) 0√ (2 − 2)
El resultado del anterior problema nos proporciona 2 formas distintas de evaluar el laplaciano, por tanto, cualquier promediado de las dos expresiones también es una discretización del laplaciano, es decir: ∆F = Fi+1,j+1 + Fi−1,j+1 + Fi−1,j+1 + Fi+1,j−1 − 4Fi,j =γ + 2h2 Fi+1,j + Fi−1,j + Fi,j+1 + Fi,j−1 − 4Fi,j +(1 − γ) + h2 +O(h) donde γ es un parámetro libre a elegir. La elección de dicho parámetro γ la haremos de forma que la discretización de ∆F respete lo máximo posible la invarianza por rotaciones de la función F (x, y). Para ello, consideremos una función tal que en un entorno de un punto (hi0 , hj0 ) tiene los siguientes valores: 1 0 0 1 0 0 1 0 0
Si calculamos ∆F en el punto central a través de la anterior fórmula obtenemos: 2 1 ∆F (hi0 , hj0 ) = γ 2 + (1 − γ) 2 2h h Ahora bien, si rotamos 45 grados, la función inicial en torno al punto (hi0 , hj0 ), obtenemos como imagen: 1 1 0 1 0 0 0 0 0
Si calculamos de nuevo ∆F en el mismo punto obtenemos: 1 2 ∆F (hi0 , hj0 ) = γ 2 + (1 − γ) 2 2h h Por lo tanto, si queremos que ambos valores de ∆F coincidan, debemos elegir γ = 2 . Hablando en términos 3 de teoría de la señal, el calculo de ∆F nos llevaría a convolucionar la imagen con la siguiente máscara: 1 h2
dan lugar a una discretización del gradiente tal que su norma euclídea es invariante por rotaciones de 45 grados. Problema 58 (2 puntos) Calcular una aproximación del gradiente de una función F (x, y) en el punto (x, y) = (0, 0) conociendo los siguientes valores: F (0, 0) = 0, F ( 1 , 0) = 1 , F (− 1 , 0) = − 1 , F (0, 1 ) = − 1 , F (0, − 1 ) = 2 2 2 2 2 2 2 1 1 1 1 1 1 1 2 , F ( 2 , 2 ) = 0, F (− 2 , − 2 ) = 0, F (− 2 , 2 ) = −1, F ( 1 , − 1 ) = 1. 2 2
Integración Numérica Métodos de Cuadratura de Gauss Sea f (x) una función deﬁnida en un intervalo [a, b], vamos a aproximar el valor de la integral de f (x) en [a, b] utilizando la evaluación de f (x) en ciertos puntos de [a, b]. Es decir, una fórmula de integración numérica se puede escribir como Z
Problema 59 (2 puntos) Aproximar el valor de la siguiente integral, utilizando las fórmulas de Legendre para n = 2 y n = 3: Z 1 ¡ 3 ¢ x − x4 dx
donde xk representa los puntos de evaluación de f (x) y wk el peso de cada punto de evaluación. Deﬁnición 2 Una fórmula de integración numérica se denomina exacta de orden M si, para cualquier polinomio P (x) de grado menor o igual que M, la fórmula es exacta. Es decir Z b N X P (x)dx = wk P (xk )
Problema 60 (2 puntos) Se consideran, para el intervalo [−1, 1], los puntos x0 = −0.5, x1 = 0 y x2 = 0.5 y los pesos w0 = w1 = w2 = 2/3. Estos puntos y estos pesos se utilizan para aproximar la integral de una función en [−1, 1]. Usar esta fórmula de integración para calcular númericamente la siguiente integral y compararla con el resultado análitico (exacto). Z
Deﬁnición 3 Se denominan polinomios de Legendre Ln (x) a la familia de polinomios dada por L0 (x) = 1, L1 (x) = x, y para n = 2, 3, .... nLn (x) = (2n − 1)xLn−1 (x) − (n − 1)Ln−2 (x)
Problema 61 (2 puntos) Encontrar, utilizando los ceros y pesos asociados a los polinomios de Legendre, cuál sería la fórmula de integración numérica de Legendre utilizando un sólo punto de interpolación. ¿Cuál sería su exactitud? Problema 62 (2 puntos) A partir de los ceros y de los pesos asociados a los polinomios de Legendre, y dado un intervalo [a, b] cualquiera, encontrar los puntos xk y los pesos wk que hacen exacta hasta orden 2N −1 una fórmula de integración numérica sobre el intervalo [a, b]. Problema 63 (2 puntos) Utilizar el resultado del problema anterior para calcular de forma exacta la siguiente integral: Z 1 ¡ 2 ¢ x − x3 dx
Teorema 14 Sean{˜k }k=1,..,N los ceros del polinomio de x Legendre LN (x). Si deﬁnimos wk = ˜ Z
entonces la fórmula de integración numérica generada por ˜ los puntos xk y los pesos wk es exacta hasta el orden 2N −1 ˜ para el intervalo [−1, 1]. Demostración [Hu] Pg. 205-209 Ejemplo 12 A continuación se exponen algunos valores de raíces xk y coeﬁcientes wk en función del grado del ˜ ˜ polinomio Ln (x) : n xk ˜ wk ˜ 2 0.5773502692 1. −0.5773502692 1 3 0.7745966692 0.5555555556 0. 0.8888888889 − 0.7745966692 0.5555555556 4 0.8611363116 0.3478548451 0.3399810436 0.6251451549 −0.3399810436 0.6251451549 − 0.8611363116 0.3478548451
Cuando el intervalo [a, b] es inﬁnito, es decir, a = −∞ o b = ∞, hay que emplear otros métodos para aproximar las integrales. En el caso [a, b] = (−∞, ∞), se utilizan los ceros de los denominados polinomios de Hermite, deﬁnidos como H0 (x) = 1, H1 (x) = 2x, y Hn (x) = 2xHn−1 (x) − 2(n − 1)Hn−2 (x) para n ≥ 2. En este caso, la fórmula de integración numérica aproxima la integral de la siguiente forma: Z
Teorema 15 Si xk son los ceros del polinomio de Hermite ˜ y deﬁnimos Z ∞ Πi6=k (x − xi) −x2 ˜ wk = ˜ e dx Πi6=k (xk − xi) ˜ −∞ entonces la fórmula de integración numérica generada por ˜ los puntos xk y los pesos wk es exacta hasta orden 2N − 1 ˜ para el intervalo (−∞, ∞). Demostración [Hu] Pg. 213-214 Ejemplo 13 A continuación se exponen algunos valores de raíces xk y coeﬁcientes wk en función del grado del ˜ ˜ polinomio Hn (x) : n xk ˜ wk ˜ 1 0. 1. 772 453 851 2 −0. 707 106 781 0. 886 226 925 5 0. 707 106 781 0. 886 226 925 5
Ejemplo 14 A continuación se exponen algunos valores de raíces xk y coeﬁcientes wk en función del grado del ˜ ˜ polinomio Ln (x) : n xk ˜ wk ˜ 1 1. 1. 2 0. 585 786 438 0. 853 553 390 3 3. 414 213 562 0. 146 446 609 3
utilizando los polinomios de Laguerre.
donde a es un número real cualquiera.
utilizando los polinomios de Hermite.
Fórmulas de Integración Numérica Compuestas Con las fórmulas que hemos visto hasta ahora, para aumentar la precisión es necesario aumentar el grado de los polinomios, lo cual resulta complejo para valores grandes de N. Una alternativa consiste en dividir previamente la integral en subintegrales de la manera siguiente: Z
Problema 65 (2 puntos) Aproximar, utilizando dos puntos de aproximación, el valor de la integral: Z ∞ 1 dx 1 + x2 −∞
Para el intervalo (0, ∞), se utilizan los polinomios de Laguerre Ln (x), deﬁnidos por L0 (x) = 1, L1 (x) = 1 − x, y Ln (x) = (2n − 1 − x)Ln−1 (x) − (n − 1)2 Ln−2 (x). para n ≥ 2. En este caso, la fórmula de integración numérica aproxima: Z
donde a = x0 < x1 < .... < xM+1 = b. A continuación se aproxima numéricamente cada una de las integrales Z xk+1 f (x)dx
Para ello, se pueden utilizar los desarrollos a partir de los polinomios de Legendre, o bien las fórmulas más simples siguientes: Fórmula del rectángulo µ ¶ Z xk+1 xk + xk+1 f (x)dx ≈ f (xk+1 − xk ) 2 xk Esta fórmula se obtiene fácilmente aproximando f (x) por el polinomio interpolador en x = xk +xk+1 . Es decir: 2 Z
Teorema 16 Si xk son los ceros del polinomio de La˜ guerre y deﬁnimos Z ∞ Πi6=k (x − xi) −x ˜ wk = ˜ e dx Πi6=k (xk − xi) ˜ 0 entonces la fórmula de integración numérica generada por ˜ los puntos xk y los pesos wk es exacta hasta orden 2N − 1 ˜ para el intervalo (0, ∞). Demostración [Hu] Pg. 211-213
Fórmula del trapecio Z xk+1 f (xk+1 ) + f (xk ) (xk+1 − xk ) f (x)dx ≈ 2 xk Esta fórmula se deduce aproximando f (x) por su polinomio interpolador en xk y xk+1 . Es decir: f (x)dx ≈ xk µ ¶ xk+1 x − xk+1 x − xk + f (xk+1 ) ≈ f (xk ) dx = xk − xk+1 xk+1 − xk xk f (xk+1 ) + f (xk ) (xk+1 − xk ) = 2 Z Z
Problema 68 (2 puntos) Aproximar, por el método de Simpson, la integral Z 1 ¡ 3 ¢ x − x4 dx
utilizando únicamente el valor de la función en los puntos: −1, − 1 , 0, 1 y 1. 2 2
Práctica 5 (Implementación Método de Integración de Simpson, 2 horas)
Ahora bien, teniendo en cuenta los resultados de la sección anterior sobre derivación numérica f 00 (xm ), se puede aproximar como f 00 (xm ) ≈ f (xk+1 ) − 2f (xm ) + f (xk ) ³ ´2
Crear una función en fortran 77 donde se implemente el método de Simpson. Los parámetros de la función serán: Los límites del intervalo de integración en precisión real y el número de subintervalos en los que se dividirá el interFórmula de Simpson valo inicial. La función a integrar se deﬁnirá aparte (como ³ ´ Z xk+1 en el caso del mérodo Müller). La función devolverá el f (xk+1 ) + f (xk ) + 4f xk +xk+1 2 f (x)dx ≈ (xk+1 − xk ) valor de la integral obtenido. Probar el método para aprox6 xk imar las siguientes integrales con diferentes valores para el parámetro de número de subintervalos y comprobar que el Esta fórmula se deduce aproximando f (x) por su deresultado se aproxima al valor exacto de la integral. sarrollo en serie de Taylor centrado en el punto xm = xk +xk+1 . Es decir: 2 Rπ 1. 0 sin(x)dx = 2 Z xk+1 f (x)dx ≈ xk R1 x ¶ Z xk+1 µ 2. 0 √1−x2 dx = 1 f 00 (xm ) 0 2 ≈ f (xm ) + f (xm )(x − xm ) + (x − xm ) dx = 2 xk µ ¶3 R∞ √ 2 f 00 (xm ) xk+1 − xk 3. −∞ e−x dx = π = 1. 772 5 = f (xm )(xk+1 − xk ) + 3 2 Nota: Las integrales con límites inﬁnitos se aproximarán cambiando el inﬁnito por un número grande. Integración numérica en dimensiones superiores En esta sección, estudiaremos las técnicas de integración numérica sobre dominios Ω de dimensión superior a 1. Para simpliﬁcar la exposición, supondremos que la dimensión es 2. Es decir, pretendemos aproximar Z F (x, y)dxdy
Por tanto, sustituyendo este valor en la aproximación anterior obtenemos Z
Aproximaremos esta integral a través de la fórmula numérica: Z X F (x, y)dxdy ≈ wij F (xi , yj )
Ω i,j
Aunque estas fórmulas sean menos precisas que las deducidas a partir de los ceros de los polinomios de Legendre, tienen la ventaja de que pueden ser utilizadas cuando sólo conocemos la función a integrar en un conjunto equiespaciado de puntos, es decir, cuando sólo conocemos f (x) en un conjunto de la forma xk = x0 + hk. Nótese que, en este caso, la integración a partir de los ceros de los polinomios de Legendre no puede utilizarse.
donde debemos elegir los puntos (xi , yj ) y los pesos wij . Para realizar esta elección se utilizan técnicas de cuadratura. Es decir, se exige que la fórmula sea exacta para polinomios en x e y de hasta un cierto grado: Z X m n xm y n dxdy = wij (xi ) (yj )
donde m y n determinan el grado de los polinomios. De estas relaciones se puede deducir, en general, los valores de los puntos y los pesos. Un caso particularmente sencillo es cuando Ω es un rectángulo [a, b]x[c, d]. En este caso, podemos escribir: Z xm y n dxdy =
AREA(T ) el área del triángulo T. En función de los vértices, el área viene determinada por ¯⎞ ⎛¯ ¯ 1 1 1 ¯ ¯ ¯ 1 AREA(T ) = ABS ⎝¯ x0 x1 x2 ¯⎠ ¯ ¯ 2 ¯ y0 y1 y2 ¯
y, por tanto, la exactitud en dimensión 2 la podemos deducir a partir de la exactitud en dimensión 1, que, en este caso, viene dada, como hemos visto anteriormente, por los polinomios de Legendre. Problema 69 (3 puntos) Deducir la fórmula de integración numérica sobre el rectángulo [−1, 1]x[−1, 1] resultante de aplicar la integración numérica en una variable en los intervalos [−1, 1], y [−1, 1]. Problema 70 (2 puntos) Deducir la fórmula de integración numérica sobre un rectángulo [a, b]x[c, d] resultante de aplicar la integración numérica en una variable en los intervalos [a, b], y [c, d].
A continuación presentaremos algunas fórmulas de integración numérica sobre triángulos utilizando diferentes números de puntos Integración sobre triángulos utilizando un punto. µ ¶ Z ∼ F x0 + x1 + x2 , y0 + y1 + y2 AREA(T ) F (x, y) = 3 3 T Integración sobre triángulos utilizando 3 puntos. Z F (x, y) ∼ AREA(T ) =
wk F (ek , yk ) x e
x2 e tos. x3 e
utilizando integración numérica. Nótese que, al igual que en dimensión 1, también podemos extender los resultados al caso en que los intervalos sean inﬁnitos, de tal forma que podemos construir fácilmente fórmulas de integración numérica para las integrales Z Z
Integración sobre triángulos utilizando 4 punZ F (x, y) ∼ AREA(T ) =
F (x, y)e−x
F (x, y)e−x−y dxdy
Problema 72 (2 puntos) Calcular una aproximación numérica de la integral Z ∞Z 2 x dxdy y2 −∞ 0 1 + e utilizando la evaluación de F (x, y) en 4 puntos. En el caso de que Ω sea un triángulo, el cálculo es un poco más complejo. Consideremos un triángulo T de vértices (x0 , y0 ), (x1 , y1 ), (x2 , y2 ). Denotaremos por
Problema 73 (2 puntos) Se considera el triángulo T de vértices (0, 0), (1, 0) y (0, 1). Deducir cual debe ser el punto (x0 , y0 ) y el peso w0 para que la fórmula de integración numérica: Z F (x, y)dxdy ≈ F (x0 , y0 )w0
sea exacta para polinomios de grado 1 en x e y. Es decir P (x, y) = ax + by + c.
Problema 74 (2 puntos) Calcular una aproximación numérica de la integral Z x2 ydxdy
Problema 77 (2 puntos) Tomar N = 2 y dibujar el lugar geométrico de los vectores x = (x1 , x2 ) que veriﬁcan que 1. k x k1 < 1 2. k x k2 < 1 3. k x k∞ < 1 Problema 78 (2 puntos) Tomar N = 2 y demostrar la siguiente desigualdad: k x k∞ ≤k x k2 ≤k x k1 Dada una matriz A de dimensión N xN , se podría deﬁnir su norma considerando la matriz como un vector de dimensión N xN . Sin embargo, resulta más útil deﬁnir la norma de una matriz subordinándola a la norma de un vector de la siguiente manera: Deﬁnición 5 Sea A una matriz y sea k . k una norma vectorial. Se deﬁne la norma de A, subordinada a la norma vectorial k . k como k A k= sup
donde Ω es el triángulo de vértices (0, 0), (2, 0) y (0, 2), utilizando 1 punto, 3 puntos y 4 puntos. ANÁLISIS NUMÉRICO MATRICIAL II En esta sección veremos algunos aspectos más avanzados del análisis matricial, incluyendo técnicas iterativas de resolución de sistemas de ecuaciones y cálculo de autovalores.
Normas de vectores y matrices Deﬁnición 4 Una norma k . k es una aplicación de un espacio vectorial E en R+ ∪ {0} que veriﬁca las siguientes propiedades: • k x k= 0 si y sólo si x = 0 • k λx k=| λ |k x k para todo λ ∈ K y x ∈ E • k x + y k≤k x k + k y k para todo x, y ∈ E. Básicamente, una norma mide la magnitud o tamaño de un vector x. Por ejemplo, en el espacio vectorial de los números reales, la norma ”natural” es el valor absoluto. Sin embargo, cuando trabajamos en varias dimensiones, esto es, x = (x1 , x2 , ...., xN ), existen múltiples formas de deﬁnir una norma. La deﬁnición más utilizada es la denominada norma p, donde p es un número real positivo, que viene deﬁnida por ÃN X
La propiedad fundamental que veriﬁca una norma matricial deﬁnida de esta forma es la siguiente: Teorema 17 Sea A una matriz y k . k una norma vectorial. Entonces, para cualquier vector x se veriﬁca que k Ax k≤k A k · k x k Demostración: Si x = 0, la desigualdad es trivial. Si x 6= 0, entonces, puesto que k x k> 0, la desigualdad anterior es equivalente a k Ax k ≤k A k kxk Ahora bien, esta desigualdad es cierta por la propia deﬁnición de k A k. Problema 79 (2 puntos) Demostrar que si A y B son dos matrices de dimensión N xN, entonces, para cualquier norma de matrices subordinada a una norma vectorial, se veriﬁca k AB k≤k A k · k B k
Un caso particularmente interesante es p = 2, que corresponde a la norma euclídea. Otro caso interesante es aquél que se produce cuando hacemos tender p hacia inﬁnito, lo que da lugar a la denominada norma inﬁnito, deﬁnida por k x k∞ = max | xi |
Problema 75 (4 puntos) Tomar N = 2 , p = 2, y demostrar que la norma k x kp veriﬁca las propiedades de la deﬁnición de norma. Problema 76 (3 puntos) Demostrar que Limp→∞ k x kp = max | xi |
A continuación veremos la relación que existe entre la norma de una matriz y sus autovalores. Empezaremos recordando algunos conceptos relacionados con los autovalores.
Deﬁnición 6 Un autovalor de A es un número λ real o complejo tal que existe un vector x, denominado autovector, tal que Ax = λx
xi , el vector x se podrá expresar como una combinación lineal de autovectores, de la forma: x = η 1 x1 + η2 x2 + .. + η N xN Al hacer Ax, y puesto que los xi son autovectores, obtenemos que
Deﬁnición 7 Se denomina polinomio característico P (λ) de la matriz A, al polinomio dado por el determinante P (λ) =| A − λI | que
Ax = η1 λ1 x1 + η2 λ2 x2 + .. + η N λN xN Como los autovectores son ortonormales, se cumple q k x k2 = (η 1 )2 + .. + (ηN )2 q 2 2 k Ax k2 = (η1 λ1 ) + .. + (ηN λN ) k Ax k2 ≤ ρ(A) k x k2 para cualquier vector x. En consecuencia, al tomar el supremo en x, la desigualdad se mantiene, lo que demuestra que k A k2 ≤ ρ(A)
Problema 80 (1 punto) Demostrar que los autovalores de A son los ceros del polinomio característico P (λ). Deﬁnición 8 Se deﬁne el radio espectral de una matriz A como ρ(A) = max{| λi | : λi autovalor de A}
Teorema 18 Sea A una matriz y k . k una norma vectorial. Entonces k A k≥ ρ(A) Demostración: Si λ es un autovalor de A, entonces existe un autovector x tal que Ax = λx, por tanto k Ax k k λx k = = |λ| ≤k A k kxk kxk Lo que demuestra el teorema. Teorema 19 Si los autovectores de una matriz A de dimensión N xN forman una base ortonormal de RN , entonces k A k2 = ρ(A) Demostración: Recordamos, en primer lugar, que una base ortonormal de vectores es un conjunto de vectores tales que cualquier otro vector se puede expresar como combinación lineal de ellos y, además, su producto escalar veriﬁca que (xi , xj ) =
Teorema 20 Si una matriz A de dimensión N xN es simétrica, entonces todos sus autovalores son reales y, además, sus autovectores forman una base ortonormal de RN . Demostración: [La-Th] Pg. 53. Problema 81 (2 puntos) Calcular los autovectores de la matriz ⎛ ⎞ 1 1 0 ⎝ 1 1 0 ⎠ 0 0 2
y determinar una base ortonormal de R3 compuesta por autovectores de A. Teorema 21 Sea A una matriz cualquiera, entonces p • k A k2 = ρ(t AA) P • k A k1 = maxj ( i | aij |) ´ ³P • k A k∞ = maxi | aij | j
donde (xi )k indica la coordenada k-ésima del vector xi . Vamos a demostrar la desigualdad k A k2 ≤ ρ(A) Dado que el teorema anterior determina la desigualdad en el otro sentido, tendríamos la igualdad, y por tanto el resultado del Teorema. Sea x un vector cualquiera. Puesto que A posee una base ortonormal de autovectores
Demostración: [La-Th] Pg. 73,75.
Problema 82 (2 puntos) Calcular las normas 2, 1 e inﬁnito de la matriz µ ¶ 1 0 A= 1 1
Problema 83 (2 puntos) Demostrar la siguiente igualdad: ρ(t AA) = ρ(A ·t A) Teorema 22 Sea A una matriz cualquiera, entonces Limn→∞ k An k= 0 ⇐⇒ ρ(A) < 1 Demostración:[La-Th] Pg. 80. Condicionamiento de una matriz El condicionamiento de una matriz es un número que nos indica la ”bondad” o buen comportamiento numérico de la matriz cuando se trabaja con ella numéricamente. Para ilustrar de qué estamos hablando, veamos el siguiente ejemplo: Ejemplo 15 Consideremos el siguiente sistema de ecuaciones ⎛ ⎞⎛ ⎞ ⎛ ⎞ 10 7 8 7 x 32 ⎜ 7 5 6 5 ⎟ ⎜ y ⎟ ⎜ 23 ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟ ⎝ 8 6 10 9 ⎠ ⎝ z ⎠ = ⎝ 33 ⎠ 7 5 9 10 v 31
Demostración: Como A(u + δu) = b + δb y Au = b, se obtiene que Aδu = δb, de donde δu = A−1 δb y, por tanto, ° ° kδuk ≤ °A−1 ° kδbk Por otro lado, también se cumple que kbk = kAuk ≤ kAk kuk de donde obtenemos que 1 kAk ≤ kuk kbk Así, multiplicando esta desigualdad con la anteriormente obtenida para kδuk , concluimos la demostración del teorema. Problema 84 (2 puntos) Demostrar que, si los autovectores de una matriz A de dimensión N xN forman una base ortonormal de RN , entonces, para la norma 2, se cumple que maxi {| λi |} χ(A) =k A k2 · k A−1 k2 = mini {| λi |} Nota: En el caso del ejemplo 15 los autovalores de la matriz son 0.01, 0.84, 3.86, y 30.29, por tanto el condicionamiento sería χ(A) = 30.29 = 3029 0.01
cuya solución es (1, 1, 1, 1). Vamos a considerar ahora el mismo sistema, perturbando ligeramente el término independiente: ⎛ 10 ⎜ 7 ⎜ ⎝ 8 7 ⎞⎛ 7 8 7 x 5 6 5 ⎟⎜ y ⎟⎜ 6 10 9 ⎠ ⎝ z 5 9 10 v ⎞ ⎞ 32.1 ⎟ ⎜ 22.9 ⎟ ⎟=⎜ ⎟ ⎠ ⎝ 33.1 ⎠ 30.9 ⎛
La solución de este sistema es (9.2, −12.6, 4.5, −1.1). Como podemos observar, a pesar de que la perturbación del sistema es del orden de 0.1, la perturbación de la solución del sistema puede llegar a ser del orden de 13.6. Consideremos de forma genérica un sistema de ecuaciones de la forma Au = b y, al mismo tiempo, el sistema de ecuaciones perturbado A (u + δu) = b + δb Nosotros queremos controlar el error relativo en la solución del sistema a partir del error relativo en el término independiente b. Es decir, queremos encontrar una estimación del tipo k δu k k δb k ≤ χ(A) kuk kbk donde χ(A) es un número que llamaremos condicionamiento de la matriz. Obviamente, cuanto más pequeño sea χ(A), mejor comportamiento numérico tendrá la matriz A.
lo cual indica un condicionamiento bastante malo. Problema 85 (2 puntos) Calcular el condicionamiento para la norma 2, de las siguientes matrices: ⎛ ⎞ 2 2 −2 A=⎝ 2 1 1 ⎠ −2 1 1 ⎛ ⎞ 2 −1 0 A = ⎝ −1 2 −1 ⎠ 0 −1 2 Cálculo de autovalores y autovectores En esta sección veremos algunos métodos elementales para el cálculo de autovalores y autovectores de matrices.
Método de Jacobi Este método se aplica a matrices reales y simétricas. Se basa en el hecho de que, dadas dos matrices A y R, se veriﬁca que los autovalores de A son los mismos que los autovalores de R−1 AR. Este método intenta diagonalizar A realizando transformaciones del tipo R−1 AR. Problema 86 (2 puntos) Sean las matrices A y R. Demostrar que la matriz A y la matriz B = R−1 AR poseen los mismos autovalores Problema 87 (2 puntos) Se considera la matriz µ ¶ 1 1 A= 1 1 calcular el ángulo α tal que la matriz µ ¶ cos α sin α R= − sin α cos α veriﬁque que la matriz B = R−1 AR sea diagonal. En el método de Jacobi se utilizan las denominadas matrices de rotación, que tienen la forma siguiente: ⎞ ⎛ 1 0 0 0 0 0 0 ⎜ 0 1 . . . . 0 ⎟ ⎟ ⎜ ⎜ 0 . cos α . sin α . 0 ⎟ ⎟ ⎜ . 1 . . 0 ⎟ Rpq (α) = ⎜ 0 . ⎟ ⎜ ⎜ 0 . − sin α . cos α . 0 ⎟ ⎟ ⎜ ⎝ 0 . . . . 1 0 ⎠ 0 0 0 0 0 0 1
Es decir, cot(2α) = cos 2α (aqq − app ) = sin 2α 2apq
Para convertir en 0 el elemento a12 = −1 de la matriz, debemos elegir α tal que cot(2α) = (a22 − a11 ) =0 2a12
t y al hacer la operación R12 AR12 obtenemos ⎛ ⎞ 1.0 0 −. 707 107 ⎝ 0 3.0 −. 707 107 ⎠ −. 707 107 −. 707 107 2.0
De donde α = − π . Por tanto, la matriz R12 es 4 ⎞ ⎛ √ 1 1 − √2 0 2 1 1 √ 0 ⎠ R12 = ⎝ √
Para evitar tener que evaluar funciones trigonométricas, que son costosas computacionalmente, y simpliﬁcar el algoritmo, podemos apoyarnos en las igualdades trigonométricas dadas en el siguiente problema: Problema 88 (3 puntos) Demostrar las siguientes igualdades trigonométricas: q tan(α) = − cot(2α) + sign(cot(2α)) 1 + cot2 (2α) ¡ ¢ donde α ∈ − π , π , sign(x) = 1 si x ≥ 0 y sign(x) = −1 4 4 si x < 0, cos α = 1 p 1 + tan2 (α) sin α = tan(α) cos α
donde los cosenos y senos están situados en las columnas y ﬁlas p y q. Al ser una matriz de rotación, se veriﬁca −1 que (Rpq (α)) =t Rpq (α). Al realizar la operación A0 =t Rpq (α)ARpq (α), sólo se ven afectadas las ﬁlas y columnas de índices p y q. Además, la matriz A0 también es simétrica. Concretamente, si A es una matriz simétrica, los cambios que se producen en A0 son los siguientes: a0 pq a0 pp a0 qq a0 pj a0 qj (app − aqq ) sin 2α + apq cos 2α 2 = app cos2 α + aqq sin2 α − apq sin 2α = = app sin2 α + aqq cos2 α + apq sin 2α = apj cos α − aqj sin α j 6= p, q = apj sin α + aqj cos α j 6= p, q
Utilizando las anteriores igualdades trigonométricas, la transformación de la matriz A mediante el método de Jacobi se puede escribir como a0 pq a0 pp a0 qq a0 pj a0 qj = = = = = 0 app − tan(α)apq aqq + tan(α)apq apj cos α − aqj sin α j 6= p, q apj sin α + aqj cos α j 6= p, q
El método de Jacobi se basa en ir modiﬁcando la matriz A mediante el procedimiento anterior, haciendo 0 los elementos no diagonales mayores en módulo. Para anular un valor a0 , basta con elegir α tal que pq (app − aqq ) sin 2α + apq cos 2α = 0 2
Problema 89 (3 puntos) Dentro del método de Jacobi para el cálculo de autovalores, demostrar las igualdades siguientes: a0 pq a0 pp a0 qq a0 pj a0 qj = = = = = 0 app − tan(α)apq aqq + tan(α)apq apj cos α − aqj sin α apj sin α + aqj cos α
j 6= p, q j 6= p, q
Veamos ahora la convergencia del método de Jacobi para el cálculo de autovalores. Teorema 24 Sea una matriz A simétrica. Sea A1 = A, y sea Ak la matriz transformada de Ak−1 , haciendo cero el elemento no diagonal mayor en módulo de la matriz Ak−1 , entonces los elementos diagonales de la matriz Ak convergen (k → ∞) hacia los autovalores de la matriz A. Además los elementos no diagonales de A convergen hacia 0. Demostración: [La-The] Pg. 576-577. Algoritmo del Método de Jacobi para el Cálculo de autovalores. Los parámetros de entrada son la matriz simétrica A, su dimensión DIM, el número máximo de iteraciones N max y la tolerancia T OL para decidir cuándo son ceros los elementos no diagonales. PARA n = 1, .., N max HACER p=2 q=1 R = ABS(A(p, q)) PARA i = 3, ..., DIM HACER PARA j = 1, ..., i − 1 HACER IF ABS(A(i, j)) > R HACER R = ABS(A(i, j)) p=j q=i FIN IF FIN PARA j FIN PARA i IF R < T OL HACER PROCEDIMIENTO TERMINADO CORRECTAMENTE. LOS AUTOVALORES SE ENCUENTRAN EN LA DIAGONAL DE A. SALIR FIN IF C = (A(q, q) − A(p, p))/(2 ∗ A(p, q)) IF C < 0 HACER T = −C − SQRT (1. + C ∗ C) ELSE T = −C + SQRT (1. + C ∗ C)
FIN IF CO = 1./SQRT (1. + T ∗ T ) SI = CO ∗ T PARA j = 1, ..., DIM HACER IF ( j 6= p AND j 6= q) HACER D = A(p, j) A(j, p) = A(p, j) = CO ∗ D − SI ∗ A(q, j) A(j, q) = A(q, j) = SI ∗ D + CO ∗ A(q, j) FIN IF FIN PARA j A(p, p) = A(p, p) − T ∗ A(p, q) A(q, q) = A(q, q) + T ∗ A(p, q) A(p, q) = A(q, p) = 0 FIN PARA n PROCEDIMIENTO TERMINADO INCORRECTAMENTE NÚMERO DE ITERACIONES MÁXIMO EXCEDIDO Veamos ahora cómo podemos calcular los autovectores. Al utilizar el método de Jacobi, vamos transformando la matriz A multiplicándola por una secuencia de matrices de rotación R1 , ...., RM , de tal forma que
−1 −1 RM · .... · R1 AR1 · .... · RM = D
donde D es una matriz diagonal que contiene los autovalores de A en la diagonal. Denotemos por B la matriz B = R1 · .... · RM . Despejando de la anterior igualdad obtenemos que AB = BD Si denotamos por bi el vector columna i de la matriz B, de la expresión anterior se obtiene que Abi = dii bi Es decir, bi es el autovector de A asociado al autovalor dii . Por tanto, la matriz B determina los autovectores. Numéricamente, para calcular la matriz B en el algoritmo anterior que calcula los autovalores, añadiremos en cada iteración las operaciones necesarios para ir obteniendo B. En primer lugar, inicializamos B a la identidad antes de entrar en el bucle. A continuación, en cada iteración haremos B = B · Ri .Ahora bien, como Ri es una matriz de rotación del tipo Rpq (α), cuando multiplicamos una matriz B por la derecha por una matriz del tipo Rpq (α) (denotemos por B 0 = B · Rpq (α) el resultado de la multiplicación) podemos observar que lo único que cambia en B son los vectores columnas p y q, que se transforman de la siguiente manera: b0 ip b0 iq = cos(α)bip − sin(α)biq = sin(α)bip + cos(α)biq i = 1, .., N i = 1, .., N
Nota. Para no tener que buscar en cada paso el máximo de los elementos no-diagonales de Ak , el algoritmo de Jacobi se puede modiﬁcar haciendo cero el primer elemento apq que se encuentre que veriﬁque |apq | ≥ T OL. Práctica 6 (Método de Jacobi para el cálculo de autovalores y autovectores 6 horas) Desarrollar en Fortran 77 la siguiente función : • FUNCTION ERROR_VECTORES(U,V,N) : Devuelve la diferencia entre los vectores U y V, de dimensión N, utilizando la fórmula :
N 1 X ABS(U (i) − V (i)) ERROR_V ECT ORES = N i=1 ABS(U (i)) + 1.
3. Las matrices de dimensión la asignatura. ⎛ 0 1 6 0 0 0 ⎜ 1 0 2 7 0 0 ⎜ ⎜ 6 2 0 3 8 0 4. A = ⎜ ⎜ 0 7 3 0 4 9 ⎜ ⎝ 0 0 8 4 0 5 0 0 0 9 5 0 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠
⎧⎛ ⎧⎛ ⎞⎫ ⎞⎫ 1 ⎬ ⎨ −1 ⎬ ⎨ √ ⎝ 0 ⎠ ↔ 2, ⎝ − 2 ⎠ ⎭ ⎩ ⎭ ⎩ 1 1 ⎧⎛ ⎞⎫ 1 ⎬ √ √ ⎨ √ 2, ⎝ 2 ⎠ ↔ 2 − 2 ⎭ ⎩ 1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎛
• FUNCTION JACOBI(A,B,N,Nmax,TOL,Niter): Realiza el cálculo de los autovalores y autovectores de una matriz simétrica A por el método de Jacobi. B es una matriz donde se guardan los autovectores por columnas. La función devuelve 0 si termina bien y −1 en caso contrario. • FUNCTION ERROR_AUTOVECTORES (A,AUTOVECTORES,AUTOVALORES,N,Nmax): Para comprobar que los autovalores λi y su autovectores xi están bien estimados, comparar ¯ para cada autovalor λi , utilizando la función ERROR_VECTORES(), los vectores A¯i y λi xi . x ¯ Devolver la expresión ERROR_AU T OV ECT ORES =
max ERROR_V ECT ORES(A¯i , λ¯i ) x x
⎞ −. 593 8 ⎜ −1. 562 8 ⎟ ⎟ ⎜ ⎜ 0.05 170 ⎟ ⎟ ⎜ 2. 11 ↔ ⎜ ⎟ ⎜ −. 400 89 ⎟ ⎠ ⎝ 1.0 . 659 86 ⎛ . 521 32 ⎜ . 741 24 ⎜ ⎜ −. 998 09 −10. 06 ↔ ⎜ ⎜ −. 855 2 ⎜ ⎝ 1.0 . 268 11
⎜ ⎜ ⎜ 16. 6 ↔ ⎜ ⎜ ⎜ ⎝ ⎛
. 363 27 . 626 28 . 901 1. 176 3 1.0 . 938 54
Comprobar los resultados obtenidos con los siguientes ejemplos, tomando T OL = 0.0001 y N iter = 1000: ⎛ ⎞ 2 2 −2 1. A = ⎝ 2 1 1 ⎠ −2 1 1 Resultado: ⎫ ⎧⎛ ⎧⎛ ⎞ ⎞⎫ 1 ⎬ ⎨ −2 ⎬ ⎨ ⎝ −1 ⎠ ↔ −2, ⎝ −1 ⎠ ↔ ⎩ ⎭ ⎩ ⎭ 1 1 ⎞⎫ ⎧⎛ ⎨ 0 ⎬ 4, ⎝ 1 ⎠ ↔ 2 ⎩ ⎭ 1 ⎛ ⎞ 2 −1 0 2. A = ⎝ −1 2 −1 ⎠ 0 −1 2 Resultado:
Nota: Obsérvese, al comparar los resultados, que los autovectores están deﬁnidos módulo la multiplicación por una constante.
⎞ 2. 984 4 ⎜ −. 796 84 ⎟ ⎟ ⎜ ⎜ 3. 088 8 ⎟ ⎟ 5. 942 ↔ ⎜ ⎜ −1. 985 4 ⎟ ⎟ ⎜ ⎠ ⎝ 1.0 −2. 165 3 ⎞ ⎛ −1. 746 5 ⎜ . 853 75 ⎟ ⎟ ⎜ ⎜ . 575 36 ⎟ ⎟ ⎜ − 2. 465 ↔ ⎜ ⎟ ⎜ −. 215 57 ⎟ ⎠ ⎝ 1.0 −1. 241 2 ⎞ ⎞ ⎛ . 729 18 ⎟ ⎜ −. 896 18 ⎟ ⎟ ⎟ ⎜ ⎟ ⎟ ⎜ ⎟ − 12. 12 ↔ ⎜ −1. 324 6 ⎟ ⎟ ⎜ 1. 826 8 ⎟ ⎟ ⎟ ⎜ ⎠ ⎠ ⎝ 1.0 −1. 767 6
Método de la potencia Teorema 25 Sea una matriz A que posee una base de autovectores tal que en módulo su autovalor máximo λmax es único. Sea un vector u1 no ortogonal al subespacio engendrado por los autovectores del autovalor λmax , entonces, si deﬁnimos la secuencia un = A se veriﬁca que Limn→∞ sign
¡¡ n n−1 ¢¢ u ,u k un k= λmax
¡ ¡¡ ¢¢¢n Limn→∞ sign un , un−1
Además, dicho autovector tiene norma 1. ¡¡ ¢¢ Teorema 26 sign un , un−1 es el signo del¢¢ producto ¡¡ escalar de un y un−1 ,¡¡ decir sign un , un−1 = 1 si es ¡ n n−1 ¢ ¢¢ ¡ ¢ u ,u ≥ 0 y sign un , un−1 = −1 si un , un−1 < 0. Demostración. En primer lugar, vamos a demostrar por inducción la siguiente igualdad: un+1 An u1 = kAn−1 u1 k
Limn→∞ kun k = ° ³ ° ´n−1 ³ ´n−1 ° ° λN °µ1 λmax x1 + ... + µN |λmax | xN ° |λmax | ° ° °= = |λmax | ° ³ ´n−2 ³ ´n−2 ° ° °µ1 λmax x1 + ... + µN |λλN | xN ° |λmax | ° ° max = |λmax | y por tanto ¡¡ n n−1 ¢¢ u ,u k un k= λmax ¡ ¡¡ ¢¢¢n Por otro lado, el término sign un , un−1 para n suﬁcientemente grande es 1n si λmax es positivo o (−1)n si λmax es negativo. Sean x1 , ...., xM los autovectores asociados a λmax , obtenemos que Limn→∞ sign ¡ ¡¡ ¢¢¢n Limn→∞ sign un , un−1 = un = k un k
Para n = 1 la igualdad se cumple por la deﬁnición de u2 . Supongamos que se cumple para n−1, y demostrémoslo para n: un =A k un k
µ1 x1 + ... + µM xM kµ1 x1 + ... + µM xM k
Con lo que queda demostrado este primer resultado. Por otro lado, como A posee una base de autovectores, que denotaremos por xi , y u1 no es ortogonal al espacio generado por los autovectores asociados a λmax , entonces u1 se puede escribir como u1 = µ1 x1 + ... + µN xN donde supondremos que x1 es un autovector asociado a λmax y que µ1 6= 0. Por la igualdad anteriormente demostrada obtenemos que
que es un autovector de λmax de norma 1. Problema 91 (3 puntos) Aplicar el método de la potencia para aproximar el autovalor máximo y el autovector asociado de las siguientes matrices, realizando 3 iteraciones en el método, hasta calcular u4 y partiendo de u1 = (1, 1). µ ¶ 2 1 A = 0 1 µ ¶ −3 0 A = 1 1
An−1 u1 µ λn−1 x1 + ... + µN λn−1 xN max N °= = ° 1 n−2 °µ1 λmax x1 + ... + µN λn−2 xN ° kAn−2 u1 k N ³ ´n−1 ³ ´n−1 Método de la potencia inversa λ x1 + ... + µN |λλN | xN µ1 |λmax | max max ° = |λmax | ° ³ ´n−2 ³ ´n−2 ° ° El método anterior también se puede utilizar para el cál°µ1 λmax x1 + ... + µN |λλN | xN ° culo del autovalor de módulo menor λ , teniendo en |λmax | ° ° max min cuenta que Cuando hacemos tender n hacia inﬁnito, todos los co1 cientes de la forma λmin = 0 max{λi autovalores de A−1 } µ ¶n λi Por tanto, si aplicamos el método anterior a A−1 , |λmax | obtenemos que la secuencia tienden hacia 0, salvo si λi = λmax . En este caso, dicho un−1 cociente es 1n , si λmax es positivo, o (−1)n , si λmax es negun = A−1 k un−1 k ativo. Por tanto, para n suﬁcientemente grande el signo de n n−1 λmax coincide con el signo del producto escalar (u , u ). veriﬁca que Además ¡¡ ¢¢ 1 Limn→∞ sign un , un−1 k un k= λmin un =
Limn→∞ sign
un es un autovector de λmin k un k En los casos prácticos, se evita calcular directamente A−1 , y se obtiene un resolviendo el sistema Aun = un−1 k un−1 k
¡¡ n n−1 ¢¢ u ,u
Para ello, calcular dos iteraciones del método de la potencia inversa partiendo de u1 = (1, 1, 1).
Métodos iterativos de resolución de sistemas lineales Estas técnicas consisten en transformar un sistema de la forma Au = b en una ecuación de punto ﬁjo de la forma u = Mu + c de tal manera que, al hacer iteraciones de la forma un = M un−1 + c se obtenga que un converge hacia u, la solución del sistema original. Ejemplo 17 Consideremos el sistema de ecuaciones 2x − y −x + 2y − z −y + 2z = 1 = 0 = 1
Problema 92 (2 puntos) Calcular el autovalor mayor µ ¶ 2 −1 y el autovector correspondiente de la matriz −1 1 utilizando el método de la potencia, realizando 2 iteraciones del método a partir de u1 = (1, 1) y tomando como norma kuk = maxi |ui |. Problema 93 (2 puntos) Utilizar el método de la potencia inversa para aproximar el autovalor menor de la matriz µ ¶ −2 1 A= 0 3 Llegar hasta u3 partiendo de u = (1, 1). Para autovalores que se encuentren entre λmin y λmax , se puede proceder de la manera siguiente: Se calcula primero una aproximación µ del autovalor λ de tal forma que µ se encuentre más cercano a λ que a cualquier otro autovalor. Por ejemplo, utilizando el método de Jacobi, si consideramos la matriz A0 = A − µI, donde µ es uno de los elementos diagonales de la matriz que resulta de aplicar el método de Jacobi, entonces se obtiene que el autovalor menor de A0 es justamente λ, y, por tanto, podemos aplicar el método de la potencia inversa anterior. Nótese que si el autovalor µ está calculado con mucha precisión, entonces el autovalor más pequeño de A0 está muy próximo a 0, y como el determinante de una matriz es el producto de sus autovalores, ello indicaría que el determinante de A0 estaría muy próximo a 0 y podemos tener problemas al resolver el sistema utilizado por ejemplo el método de GAUSS a través de la función de la librería an.h IGAUSS(). Para evitar esto, podemos perturbar ligeramente el valor de µ para que IGAUSS() no dé problemas. Algorítmicamente, quedaría como sigue: Si µ es el autovalor que estamos tratando, haremos 1 = 10−11 A = A − µId J = IGAU SS(A0 , ......) IF (J.N E.0) T HEN = ∗ 10. µ = µ(1 + ) GOT O 1 EN DIF
Buscar la solución de este sistema es equivalente a buscar un vector u = (x, y, z) que veriﬁque que x = y z = = 1+y 2 x+z 2 1+y 2
Hacer iteraciones de esta ecuación de punto ﬁjo consiste en partir de una aproximación inicial (x1 , y1 , z1 ) y hacer iteraciones de la forma xn yn zn = = = 1 + yn−1 2 xn−1 + zn−1 2 1 + yn−1 2
En este caso, la solución exacta del sistema es u = (1, 1, 1). Si hacemos iteraciones del esquema anterior a partir de la aproximación inicial u1 = (0, 0, 0), obtenemos que x2 y2 z2 = = = 1+0 1 = 2 2 0+0 =0 2 1+0 1 = 2 2
Como puede observarse, las sucesivas iteraciones se van aproximando a la solución u = (1, 1, 1). En este caso, la matriz M y el vector c que determinan el esquema iterativo vienen dados por ⎛ ⎛ 1 ⎞ ⎞ 0 1 0 2 2 MJ = ⎝ 1 0 1 ⎠ cJ = ⎝ 0 ⎠ 2 2 1 0 1 0 2 2 Teorema 27 Si el esquema iterativo un = M un−1 + c converge hacia un vector u, entonces u veriﬁca que u = Mu + c
De la misma forma, obtenemos que ⎛ ⎞ ⎛ ⎞ ⎛ ⎞ 0.5 .84 .98 u3 = ⎝ 0.25 ⎠ .... u8 = ⎝ .73 ⎠ .... u17 = ⎝ .96 ⎠ 0.5 .84 .98
devuelve el signo del producto escalar de los vectores uf y vf de dimensión Nf (12 líneas de código como máximo), y la función AUTOVALOR_MAXIMO(Af,uf,Nf,Nfmax,Nﬁter,Tolf ) que devuelve el autovalor máximo de una matriz y su autovector por el método de la potencia. Los parámetros son la matriz Af, el vector candidato inicial uf, Nf la dimensión real, Nfmax, la dimensión para coger memoria, Nﬁter número máximo de iteraciones, y Tolf la tolerancia. Esta función devuelve el valor 2.**120 si P termina no correctamente. Tomar como norma kuk = i ABS(ui ) (28 líneas de código como máximo). Problema 96 (2 puntos) Calcular 3 iteraciones del método de Jacobi para resolver el sistema ⎛ ⎞⎛ ⎞ ⎛ ⎞ 1 −1 0 x −1 ⎝ −1 2 0 ⎠ ⎝ y ⎠ = ⎝ 3 ⎠ 0 −1 3 z 1 partiendo de u1 = (0, 0, 0)
Método de Gauss-Seidel Existen diferentes métodos para convertir un sistema de la forma Au = b en una ecuación de punto ﬁjo u = M u + c. Todas se basan en descomponer A de la forma A = L + D + U, donde D es la matriz diagonal que corresponde a la parte diagonal de A, L es la matriz triangular inferior que corresponde a la parte de A situada por debajo de la diagonal, y U es la matriz triangular superior que corresponde a la parte de A situada por encima de la diagonal. Este método consiste en tomar MGS cGS = (D + L) = (D + L)
A efectos prácticos, la aplicación de este método no requiere el cálculo directo de la matriz inversa (D + L)−1 , puesto que el paso de una iteración a otra puede hacerse de la siguiente forma: un 1 −a12 un−1 − ... − a1N un−1 + b1 2 N a11 −a21 un − a23 un−1 ... − a2N un−1 + b2 1 3 N = a22 . −aN1 un − aN 2 un ... − aNN −1 un −1 + bN 1 2 N = aN N =
Es el que se ha utilizado en el ejemplo anterior. El paso de una iteración a otra del método de Jacobi puede expresarse de la siguiente forma: un 1 un 2 un N −a12 un−1 − ... − a1N un−1 + b1 2 N a11 −a21 un−1 − a23 un−1 ... − a2N un−1 + b2 1 3 N = a22 . −aN1 un−1 − aN 2 un−1 ... − aNN −1 un−1 + bN 1 2 N −1 = aN N =
Si hacemos un barrido para el cálculo de la solución de arriba hacia abajo, y vamos actualizando las componentes del vector aproximación según las vamos calculando, obtenemos el método de Gauss-Seidel. Por tanto, básicamente, podemos decir que la diferencia entre el método de Gauss-Seidel y el método de Jacobi es que en el método de Gauss-Seidel se actualiza el vector aproximación después del cálculo de cada componente, y en el caso de Jacobi se actualiza sólo al ﬁnal, después de haber calculado todas las componentes por separado. Ejemplo 18 Vamos a aplicar el método de Gauss-Seidel al sistema del ejemplo anterior, es decir 2x − y −x + 2y − z −y + 2z
Problema 95 (3 puntos) en Fortran las funciones SIGNO_PRODUCTO_ESCALAR(uf,vf,Nf )
Las iteraciones del método de Gauss-Seidel aplicado a este sistema consisten en xn yn zn = = = 1 + yn−1 2 xn + zn−1 2 1 + yn 2
estado de la solución en la etapa anterior, de la forma siguiente: un 1 un 2 un N −a12 un−1 − ... − a1N un−1 + b1 2 N + (1 − w)un−1 1 a11 n−1 n −a21 u1 ... − a2N uN + b2 = w + (1 − w)un−1 2 a22 . −aN1 un ... − aN N−1 un + bN 1 N−1 = w + (1 − w)un−1 N aNN = w
Si hacemos iteraciones del esquema anterior a partir de la aproximación inicial u1 = (0, 0, 0), obtenemos que x2 y2 z2 = = = 1+0 1 = 2 2 1 1 2 +0 = 2 4 1+ 1 5 4 = 2 8
La elección del parámetro w es, en general, un problema difícil. Sin embargo, en el caso de matrices tridiagonales, es decir, matrices con todos los elementos nulos salvo la diagonal principal y sus codiagonales, el siguiente resultado muestra la forma de calcular el valor óptimo de w. Teorema 28 Si A es una matriz tridiagonal y ρ(MJ ) < 1, entonces el valor de w que optimiza la velocidad de convergencia del método es: wopt = 2 p 1 + 1 − ρ(MJ )2
De la misma forma, obtenemos que ⎛ ⎞ ⎛ ⎞ .5 . 976 56 u3 = ⎝ . 25 ⎠ .... u8 = ⎝ . 976 56 ⎠ . 625 . 988 28
Como puede observarse de la expresión anterior, el valor de wopt se encuentra siempre entre 1 y 2. Demostración [La-Th]. Pg.358-362. Ejemplo 19 Vamos aplicar el método de relajación al sistema del ejemplo anterior, es decir 2x − y −x + 2y − z −y + 2z = 1 = 0 = 1
Problema 98 (2 puntos) Calcular 3 iteraciones del método de Gauss-Seidel para resolver el sistema ⎛ ⎞⎛ ⎞ ⎛ ⎞ 1 −1 0 x −1 ⎝ −1 2 0 ⎠ ⎝ y ⎠ = ⎝ 3 ⎠ 0 −1 3 z 1 partiendo de u1 = (0, 0, 0)
Problema 99 (1 punto) Una variante del método de Gauss-Seidel consiste en tomar M = (D + U )−1 (−L), y −1 c = (D + U ) b. Indicar, en este caso, qué diferencias de implementación habría con respecto al caso anterior.
1 En este caso, ρ(MJ ) = √2 y wopt = 1. 17. Las iteraciones del método de relajación aplicado a este sistema consisten en
Método de relajación El objetivo de este método es intentar mejorar el método de Gauss-Seidel introduciendo un parámetro de relajación w. Se toman, en este caso, Mw cw = (D + wL)
Si hacemos iteraciones del esquema anterior a partir de la aproximación inicial u1 = (0, 0, 0) y tomando w = wopt = 1.17, obtenemos que ⎞ ⎛ ⎞ ⎛ ⎞ . 585 . 686 . 999 u2 = ⎝ . 342 ⎠ .... u3 = ⎝ . 802 ⎠ ... u8 = ⎝ . 999 ⎠ . 785 . 921 . 999 ⎛
partiendo de u1 = (0, 0, 0). Calcular previamente el parámetro de relajación óptimo. Convergencia de los métodos iterativos Vamos a denotar por en = un − u el error relativo entre la solución del sistema u y la aproximación en la etapa n, un . Teorema 29 Se considera el esquema iterativo un = M un−1 + c. Entonces en = M n−1 e1 Demostración: La solución del sistema satisface que u = M u + c. Restando esta igualdad de la igualdad un = M un−1 + c, obtenemos que un − u = M (un−1 − u) = M n−1 (u1 − u) Teorema 30 El método iterativo un = M un−1 + c converge para cualquier aproximación inicial si y sólo si ρ(M ) < 1. Demostración: El resultado es inmediato a partir del hecho de que una matriz M n converge hacia 0 cuando n → ∞ si y sólo si ρ(M ) < 1 Teorema 31 Si en el método de relajación w ∈ (0, 2), / entonces ρ(Mw ) ≥ 1. Demostración: En primer lugar, observamos que las matrices D +Lw y (1−w)D −wU son matrices triangulares y, por tanto, su determinante es el producto de los elementos diagonales. Además, teniendo en cuenta que el determinante del producto de dos matrices es el producto de sus determinantes y que el determinante de la matriz inversa es el inverso del determinante, obtenemos que
entonces el método de Jacobi asociado al sistema Au = b converge para cualquier aproximación inicial. Demostración: En primer lugar, observamos que la matriz MJ puede expresarse como: ⎛ 0 − a12 − a13 . − a1N a11 a11 a11 a21 a23 ⎜ − a22 0 − a22 . − a2N ⎜ a22 ⎜ . . . . . ⎜ aN −1,2 a ⎜ − aN −1,1 − aN−1,N−1 . 0 − aNN−1,N ⎝ aN−1,N−1 −1,N −1 aN,1 aN,2 a − aN,N − aN,N . − N,N −1 0 aN,N
Teniendo en cuenta que las normas 1 e inﬁnito de una matriz son el máximo de las sumas por ﬁlas o columnas en valor absoluto, se tiene, por las condiciones del teorema, que kMJ k < 1 para la norma 1 o inﬁnito. Por tanto, el teorema se concluye teniendo en cuenta que cualquier norma de una matriz es siempre mayor o igual que su radio espectral. Este resultado se puede generalizar un poco al caso de matrices irreducibles de la siguiente forma: Deﬁnición 9 Una matriz A es irreducible si un sistema de la forma Au = b no puede descomponerse en dos subsistemas independientes de dimensión menor Dicho de otra forma, una matriz es irreducible si el cambio de cualquier valor del vector b del sistema Au = b afecta a todos los elementos del vector u. Teorema 33 Si A es una matriz irreducible y se veriﬁca que X | aii |≥ | aij | ∀i.
con la desigualdad estricta en al menos una ﬁla o columna, entonces los métodos iterativos convergen. Demostración. [La-The] Pg.346-347. Ejemplo 20 La matriz del sistema ejemplo tratado anteriormente, esto es ⎛ ⎞ 2 −1 0 ⎝ −1 2 −1 ⎠ 0 −1 2 satisface las hipótesis del Teorema anterior.
Por lo tanto, como el determinante de una matriz es el producto de sus autovalores, obtenemos que, si w ∈ / (0, 2), entonces |1 − w| ≥ 1 y, en consecuencia, Mw posee al menos un autovalor de módulo mayor o igual que uno. Teorema 32 Si una matriz A veriﬁca que X | aij | ∀i. | aii |>
Problema 101 (2 puntos) Escribir en Fortran la función siguiente: CONDICIONAMIENTO(Af,Nf,Nfmax,TOLf,Nﬁter) que devuelve el condicionamiento de una matriz utilizando el método de Jacobi para calcular los autovalores. Se supondrá implementada la función JACOBI(A,N,Nmax,TOL,Niter) que devuelve 0 si termina bien y 1 si termina mal. La función CONDICIONAMIENTO devuelve 2.*120 si termina mal porque Jacobi da un error o se produce una división por cero. Los parámetros son la matriz Af, Nf la dimensión real, Nfmax, la dimensión para coger memoria, Nﬁter número máximo de iteraciones, y Tolf la tolerancia (21 líneas de instrucciones como máximo). Problema 102 (2 puntos) Demostrar que, si una matriz A veriﬁca que por ﬁlas o columnas su suma es siempre igual a 0, entonces el determinante de A es cero, y por tanto el sistema asociado a A no tiene solución.
4. Los sistemas ejemplos del directorio de la asignatura. Estos ejemplos tienen siempre como solución el vector (1, 1, ...., 1).
Método de Newton-Raphson para sistemas de ecuaciones no lineales En las aplicaciones reales, muchas veces nos encontramos con sistemas no lineales de ecuaciones. Por ejemplo, calcular las raíces, reales o complejas, de un polinomio de grado 2 dado por P2 (z) = az 2 + bz + c, donde z = x + yi, es equivalente a resolver el sistema ax2 + bx − ay 2 + c = 0 2ayx + by = 0 que es un sistema no lineal de ecuaciones. En general, un sistema no lineal de ecuaciones de dimensión N, se escribe como N ecuaciones del tipo f1 (u1 , ...., uN ) = 0 f2 (u1 , ...., uN ) = 0 . fN (u1 , ...., uN ) = 0 donde f (u) = (f1 (u), f2 (u), ..., fN (u)) es una función de <N → <N , y u = (u1 , ...., uN ). El método de NewtonRaphson para sistemas de ecuaciones se basa en desarrollar por Taylor la función f y truncar el desarrollo para que quede un sistema lineal, es decir ¡ ¢ ¡ ¢ f (u) = f (u0 ) + ∇f (u0 ) u − u0 + O k u − u0 k2
Problema 103 (3 puntos) Dado un sistema iterativo un = M un−1 + c Demostrar que, aunque el radio espectral de M sea mayor que 1, si u1 y c son combinaciones lineales de autovectores de M correspondientes a autovalores de módulo menor que 1, entonces el método converge.
Práctica 7 (Método de relajación, 2 horas) Desarrollar una función en Fortran 77 donde se implemente el método de relajación. Los parámetros de la función serán: la matriz A, el vector b, un vector u donde se almacenará la solución, y que inicialmente será el vector aproximación inicial, que por defecto se tomará 0, el parámetro de relajación w, el número máximo de iteraciones N max, y la tolerancia T OL para evaluar la diferencia entre un y un−1 . La función devolverá el número de iteraciones necesarias para alcanzar la solución. Si el método no converge devuelve −1. Comparar la diferencia en la velocidad de convergencia entre el método de Gauss-Seidel y el Método de relajación. Probar el método para los sistemas 1 −1 0 x −1 1. ⎝ −1 2 0 ⎠ ⎝ y ⎠ = ⎝ 3 ⎠ 0 −1 3 z 1 ⎛ ⎛ ⎞⎛ ⎞ ⎛ ⎞
donde u0 es una aproximación de la solución de f (u) = 0. Si truncamos el desarrollo e igualamos a 0 (para aproximar la raíz) obtenemos que la raíz del sistema lineal se obtiene resolviendo el sistema ∇f (u0 )z u1 = −f (u0 ) = u0 + z
En el caso general, a partir de una aproximación un se obtiene la aproximación un+1 en dos etapas: ∇f (un )z un+1 = −f (un ) = un + z
⎞⎛ ⎞ ⎛ ⎞ 2 −1 0 x 1 2. ⎝ −1 2 −1 ⎠ ⎝ y ⎠ = ⎝ 0 ⎠ 0 −1 2 z 1 ⎞⎛ ⎞ ⎛ ⎞ 1 3 3 x 7 3. ⎝ 3 1 3 ⎠ ⎝ y ⎠ = ⎝ 7 ⎠ 3 3 1 z 7 ⎛
Si partimos de u1 = (1, 1), para obtener u2 tenemos que resolver µ ¶µ ¶ µ ¶ 2 −2 z1 −1 = 2 2 −2 z2 ¡ 3 1¢ que tiene por solución − 4 , − 4 . Por tanto, u2 viene dado por µ ¶ µ 3 ¶ µ 1 ¶ 1 −4 2 4 u = + = 3 1 −1 4 4 Para calcular u3 , tenemos que resolver el sistema µ 1 ¶µ ¶ µ ¡ ¢2 ¡ ¢2 ¶ 1 z1 −3 − 3 +1 2 2 4 4 =− 3 1 6 z2 2 2 16 ¡ ¢ 9 cuya solución es − 13 , − 40 . Por tanto, u3 viene dado por 40 µ 1 ¶ µ 13 ¶ µ ¶ 3 − 40 − 40 3 4 u = + = 3 9 39
Tomemos como aproximación inicial u1 = (1, 1). El sistema que hay que resolver para pasar de una iteración a otra es ¶µ ¶ µ 2 ¶ µ 2 z1 xn − yn + 1 2xn −2yn =− 2yn xn 2yn 2xn z2
INTERPOLACIÓN DE FUNCIONES II Esta sección es la continuación natural del tema interpolación de funciones visto anteriormente. Por motivos de coordinación entre los programas teórico y práctico de la asignatura, el tema de interpolación de funciones se dividió en dos partes, siendo ésta la segunda parte.
Interpolación de Hermite En ocasiones, resulta de interés interpolar no sólo el valor de la función en ciertos puntos {xi }i=0,..,N , sino también el valor de sus derivadas. Un ejemplo clásico de ello es el desarrollo de Taylor de una función en un punto a. En este caso, aproximamos f (x) por un polinomio de grado N , PN (x) tal que f (x) y PN (x) poseen las mismas derivadas en el punto a desde el orden 0 hasta el orden N. f 0 (a) f N ) (a) (x − a) + .. + (x − a)N 1! N!
que ya es una buena aproximación de la solución exacta dada por el vector (0, 1). Problema 104 (3 puntos) Calcular 2 iteraciones del método de Newton-Raphson no lineal para aproximar una raíz del sistema de ecuaciones x2 + y 2 − 1 = 0 y−x = 0 partiendo de (x, y) = (1, 1). Problema 105 (2 puntos) Plantear el algoritmo necesario para calcular, utilizando el método de NewtonRaphson, las raíces complejas o reales de un polinomio de grado 3.
donde ξ es un valor intermedio entre x y a. En el caso general, donde buscamos un polinomio P (x) tal que él y todas sus derivadas hasta un cierto orden M coincidan con una función f (x) en los puntos {xi }i=0,..,N , se utilizan los denominados polinomios base de Hermite Hi,j (x), que son polinomios de grado menor o igual que (N + 1)(M + 1) − 1 dados por las siguientes condiciones: ½ ∂ l Hi,j 1 si l = j y k = i (xk ) = 0 l 6= j o k 6= i ∂xl A partir de los polinomios base de Hermite, el polinomio interpolador de Hermite se deﬁne como: P (x) =
(xi )Hi,j (x)
A partir de u1 = (1, 1), calcular u2 y u3 utilizando el método de Newton-Raphson para aproximar un cero del sistema no lineal. Problema 107 (2 puntos) Calcular una iteración del método de Newton-Raphson no lineal para aproximar una raíz del sistema de ecuaciones exyz − 1 = 0 y2 − z 3 − 2 = 0 (z − 1)x4 − 3 = 0 partiendo de (x, y, z) = (1, 1, 1).
Problema 108 (3 puntos) Calcular los polinomios base de Hermite que corresponden a tomar como puntos de interpolación x0 = −1, x1 = 1, y el orden de derivación M = 1. Interpolación por splines cúbicos Uno de los problemas básicos del polinomio interpolador de Lagrange, es que, para valores grandes de N, los polinomios de grado N pueden tener un carácter fuertemente oscilante, y los resultados obtenidos por la interpolación pueden no ser muy satisfactorios, como indica el ejemplo siguiente.
Ejemplo 22 El polinomio base de Lagrange centrado en 0 sobre los puntos xi = −5, −4, −3, −2, −1, 0, 1, 2, 3, 4, 5 es (x2 − 1)(x2 − 4)(x2 − 9)(x2 − 16)(x2 − 25) P 0 (x) = −14400 Tiene un marcado carácter oscilante como muestra su gráﬁca en el intervalo [−5, 5]. Ejemplo 23 (x2 − 1)(x2 − 4)(x2 − 9)(x2 − 16)(x2 − 25) = −14400
las condiciones anteriores, entonces ai di bi = f (xi ) i = 0, .., N ci+1 − ci = i = 0, .., N − 1 3hi ai+1 − ai hi (2ci + ci+1 ) = − hi 3
i = 0, .., N − 1
hi−1 ci−1 + 2(hi−1 + hi )ci + hi ci+1 = 3(ai+1 − ai ) 3 (ai − ai−1 ) − = hi hi−1 para i = 1, .., N − 1.
i Demostración De la condición P3 (xi ) = f (xi ), se obtiene de forma inmediata que ai = f (xi ). De la condición
i ∂ 2 P3 ∂x2 (xi+1 ),
2ci+1 = 6di hi + 2ci de donde, despejando, obtenemos que
0 2.5 x 5
1.25 0 -5 -2.5
i De la Condición P3 (xi+1 ) = f (xi+1 ), se obtiene que
di h3 + ci h2 + bi hi + ai = ai+1 i i Para evitar este problema de oscilaciones de los polinomios de Lagrange, cuando se trabaja con muchos puntos de interpolación, se suele interpolar la función utilizando polinomios a trozos, deﬁniendo un polinomio distinto para cada intervalo [xi , xi+1 ]. La técnica más conocida es la interpolación por splines cúbicos, que son polinomios de grado 3. Por tanto, tendremos un polinomio de grado 3 disi tinto P3 (x) = di (x− xi )3 + ci (x−xi )2 + bi (x− xi ) +ai para cada intervalo [xi , xi+1 ]. Si hay N +1 puntos, el número de polinomios es N. Para deﬁnir estos polinomios, se imponen las siguientes condiciones:
i P3 (xi ) = f (xi ) i = 0, .., N − 1 i P3 (xi+1 ) = f (xi+1 ) i = 0, ..., N − 1
Despejando, obtenemos que bi = ai+1 − ai − di h2 − ci hi = i hi ai+1 − ai hi (2ci + ci+1 ) − hi 3
Finalmente, de la condición obtiene que
i−1 ∂P3 ∂x (xi ),
bi = 3di−1 h2 + 2ci−1 hi−1 + bi−1 i−1 y, despejando todo en función de ci , se obtiene la relación hi−1 ci−1 + 2(hi−1 + hi )ci + hi ci+1 = = 3(ai+1 − ai ) 3 (ai − ai−1 ) − hi hi−1
i+1 ∂P3 (xi+1 ) i = 0, .., N − 2 ∂x i+1 ∂ 2 P3 (xi+1 ) i = 0, ..., N − 2 ∂x2
= Vamos a introducir la notación hi xi+1 − xi . Nótese que, para deﬁnir los polinomios, tenemos que buscar 4N valores, es decir: a0 , ...., aN −1 , b0 , ..., bN−1 , c0 , ....., cN−1, d0 , ...., dN−1. Por razones técnicas, como veremos posteriormente, vamos a utilizar también los valores aN y cN .
i Teorema 34 Si una familia de polinomios P3 (x) = di (x− 3 2 xi ) + ci (x − xi ) + bi (x − xi ) + ai , i = 0, .., N, satisface
Nótese que esta última relación determina un sistema de ecuaciones donde las incógnitas son las variables ci . Dicho sistema tiene N +1 incognitas (c0 , ..., cN ) y N −1 ecuaciones. Para completar dicho sistema, hay que añadir una ecuación que involucre a c0 y otra ecuación que involucre a cN . Para añadir estas dos ecuaciones hay dos procedimientos estándares. El primero consiste simplemente en ﬁjar c0 = cN = 0, lo que signiﬁca que c0 cN
El segundo procedimiento se utiliza cuando utilizamos los valores de f 0 (a) y f 0 (b). En este caso, imponemos que
−2.2 = 1.733 3 −4.4 + 2.8 = −1 − = −0.467 3 5.6 + 0 = 2− = 0.133 3 = 1−
Por tanto, los polinomios son P1 (x) = 1.667 (x − 1)3 − 2.2 (x − 1)2 −0.467 (x − 1) + 1
P0 (x) = −0.733x3 + 1.733x
P2 (x) = −0.933 (x − 2) + 2.8 (x − 2) + 0.133 (x − 2)
Por lo tanto, siguiendo con el resultado del teorema anterior, para calcular los splines cúbicos es necesario, en primer lugar, tomar ai = f (xi ). A continuación, se resuelve un sistema de ecuaciones tridiagonal para el cálculo de los ci . Los bj y dj se calculan directamente a partir de las relaciones mostradas en el teorema anterior. Ejemplo 24 Vamos a calcular los polinomios interpoladores utilizando splines cúbicos al interpolar la función f (x) en los puntos x = 0, 1, 2, y 3, sabiendo que f (0) = 0, f (1) = 1, f (2) = 0, f (3) = 2, tomando c0 = c3 = 0. En este caso hi = 1. Debemos deﬁnir 3 polinomios distintos que corresponden a los intervalos [0, 1], [1, 2], y [2, 3]. Los términos ai vienen dados por ⎞ ⎛ 0 a0 ⎜ a1 ⎟ ⎜ 1 ⎟ ⎜ ⎜ ⎝ a2 ⎠ = ⎝ 0 2 a3 ⎛ ⎞ ⎟ ⎟ ⎠
A continuación se muestra una gráﬁca con los 3 polinomios concatenados en el intervalo [0, 3] : Como puede observarse, por las condiciones sobre las derivadas que hemos impuesto, no es posible distinguir geométricamente, al trazar la curva, cuales son los puntos de unión entre los tres polinomios. Es decir, parece, a simple vista, el trazado de una única función. Veamos ahora gráﬁcamente el perﬁl de la derivada de los polinomios P0 (x), P1 (x), y P2 (x).
2.75 2.5 2.25 2 1.75 1.5 1.25 1 0.75 0.5 0.25 0 -0.25 -0.5 -0.75 -1 -1.25
2.75 3 x
como puede observarse, tampoco sobre la derivada se aprecian los puntos de unión de los polinomios. Sin embargo, sobre la gráﬁca de la derivada segunda los puntos de unión se detectan en los lugares donde encontramos un pico, tal y como se muestra en la gráﬁca de la derivada segunda siguiente:
cuya solución es ¶ µ ¶ µ c1 −2.2 = . 2.8 c2 Los valores bi y di vienen dados por d0 d1 d2 = = = −2.2 − 0 = −0.733 3 2.8 + 2.2 = 1.667 3 0. − 2.8 = −0.933 3
y 5 4 3 2 1 0 0 -1 x -2 -3 -4 0.25 0.5 0.75 1 1.25 1.5 1.75 2 2.25 2.5 2.75 3
Problema 109 (3 puntos) Calcular los polinomios que determinan la interpolación por splines cúbicos de la fun¡ ¢ ción f (x) = sin π x para los puntos x = −1, 0, 1 y 2. 2 La interpolación a través de la función seno cardinal Una base de funciones interpolantes muy utilizada en la teoría de Fourier es la base formada a partir de la función seno cardinal, deﬁnida por sin c(x) = cuya gráﬁca es sin(x) x
0 0 0.5 1 1.5 2 2.5 x 3
Ejemplo 27 Vamos a comparar gráﬁcamente el resultado de interpolar la función del ejemplo anterior utilizando la interpolación de Lagrange normal, la interpolación por splines cúbicos y la interpolación a través de la función seno cardinal. El polinomio interpolador de Lagrange se puede calcular fácilmente y da como resultado 5 P (x) = x − x(x − 1) + x(x − 1)(x − 2) 6 En la siguiente ﬁgura se muestran juntas las gráﬁcas del polinomio de Lagrange (línea a trozos), los polinomios de la interpolación por splines cúbicos (línea sólida), y la interpolación utilizando la función sin c(x) (línea a trozos).
0.25 0 -50 -25 0 25 x 50
Esta función tiene la propiedad de que en x sin c(0) = 1, y para cualquier entero i distinto sin c(πi) = 0. Dada una función f (x), su función polante en los puntos xi = a · i para i = M, ..., N dada por la función ¢ ¡ sin(π x − i ) e ¡x a ¢ f (xi ) f (x) = π a −i i=M
= 0, de 0, interviene
1.4 1.2 1 0.8 0.6 0.4 0.2 0 0 0.25 0.5 0.75 1 1.25 1.5 1.75 2 2.25 2.5 2.75 3 x
Ejemplo 26 Consideremos la función f (x), deﬁnida en los puntos x = 0, 1, 2, y 3, tal que f (0) = 0, f (1) = 1, f (2) = 0, f (3) = 2. La interpolación de esta función utilizando la función seno cardinal viene dada por la función sin(π (x − 1)) sin(π (x − 3)) e f (x) = +2 π(x − 1) π(x − 3)
Como puede observarse, la interpolación por splines cúbicos es la menos oscilante. Por otro lado, cuando el número de puntos de interpolación aumenta, la diferencia entre los diferentes tipos de interpolación también lo hace. Problema 110 (2 puntos) Calcular la función que interpola, utilizando la función sin c(x), la función f (x) = sin(x) en los puntos x = −π, − π , 0, π , π. 2 2 La interpolación trigonométricos a través de polinomios
Dada una función f (x), deﬁnida en el intervalo [−π, π], pretendemos aproximar f (x) como f (x) ≈
donde ck son coeﬁcientes, en general complejos. El siguiente resultado determina la forma de calcular dichos coeﬁcientes ck : Teorema 35 Los coeﬁcientes ck que minimizan el error cuadrático medio !2 Z π Ã N X ikx f (x) − ck e dx E(c−N , ..., cN ) =
0.25 0 -2.5 -1.25 0 1.25 2.5 x
Demostración En primer lugar, observamos que, dada la forma cuadrática del funcional E(c−N , ..., cN ), éste debe poseer mínimos. Por otro lado, en un mínimo, las derivadas parciales de E(c−N , ..., cN ) con respecto a cualquier ck son cero, y por tanto ! Z π Ã N X ∂E ilx (c−N , ..., cN ) = cl e f (x) − eikx dx = 0 ∂ck −π
Problema 111 (3 puntos) Calcular el polinomio trigonométrico, tomando N = 2, que interpola la función f (x) = |x| en el intervalo [−π, π]. Aproximación por mínimos cuadrados La aproximación mínimo cuadrática aproxima, a través de una función, un conjunto de valores de forma global, sin exigir que la función aproximante pase exactamente por ese conjunto de puntos. Dado un conjunto de valores {(xi , yi )}i=1,..,N , la aproximación mínimo cuadrática lineal consiste en buscar la recta y = ax + b, tal que la función de error cuadrático E(a, b) = sea mínima. Teorema 36 Los valores a y b que minimizan el error cuadrático anterior son P P P N N xi yi − N xi N yi i=1 i=1 i=1 a = ´2 PN 2 ³PN N i=1 xi − i=1 xi PN PN PN 2 PN i=1 xi i=1 yi − i=1 xi yi i=1 xi b = ´2 PN 2 ³PN N i=1 xi − i=1 xi
con lo que el resultado del teorema sale de forma inmediata, teniendo en cuenta que ½ Z π 2π si l = −k eilx eikx dx = 0 si l 6= k −π
Ejemplo 28 Consideremos la función ½ 1 si x ∈ [− π , π ] 2 2 f (x) = 0 si x ∈ [− π , π ] / 2 2
Vamos a calcular el polinomio trigonométrico interpolante para N = 3. Los valores de ck son Rπ f (x)dx 1 −π c0 = = 2πR 2 π f (x)eix dx 1 c1 = c−1 = −π = π R π 2π 2ix f (x)e dx c2 = c−2 = −π =0 2π Rπ f (x)e3ix dx 1 =− c3 = c−3 = −π 2π 3π Por tanto, el polinomio trigonométrico interpolador es
Demostración En primer lugar, observamos que, dada la forma cuadrática que tiene el funcional, debe poseer un
mínimo. Además, en un mínimo del funcional E(a, b), las derivadas parciales son cero, y por tanto X ∂E (axi + b − yi ) xi = 0 (a, b) = 2 ∂a i=1
X ∂E (axi + b − yi ) = 0 (a, b) = 2 ∂b i=1
Esto da lugar a un sistema lineal de ecuaciones cuyas incógnitas son a y b, y cuya resolución lleva al resultado establecido en el teorema. Problema 112 (2 puntos) Calcular la aproximación mínimo cuadrática lineal de la tabla xi 0 1 2 3 yi 0 1 0 2
[Is-Ke] Isaacson E., Keller H. ”Analysis of Numerical Methods”. John Wiley and Sons, 1966. Uno de los libros clásicos más conocidos en Análisis Numérico. Destaca por el rigor matemático en su exposición. [Ki-Ch] Kincaid D., Cheney W. "Análisis Numérico". Addison-Wesley Iberoamericana, 1994. Excelente libro de base para un curso de Métodos Numéricos. Contiene todos los tópicos habituales con una descripción muy completa y detallada. Los algoritmos están muy bien descritos a través de un seudocódigo. Trae una buena selección de problemas. [La-Th] Lascaux P., Théodor R. "Analyse numérique matricielle appliquée à l’art de l’ingénieur. Vol. 1 Méthodes directes y Vol. 2 Méthodes itératives ", Masson, 1993. Esta obra, dividida en dos volúmenes, trata en profundidad todos los tópicos relacionados con el Análisis Numérico Matricial. Su mayor virtud es el rigor matemático con el que se tratan los temas y una cuidada presentación. [St] Stewart G.W. ”Afternotes on Numerical Analysis” SIAM, 1996. Esta obra, sin pretender ser exhaustiva, muestra las últimas tendencias en cuanto a la enseñanza de los conceptos básicos del Análisis Numérico.
BIBLIOGRAFÍA BÁSICA [Bu-Fa] Burden R., Faires D. ”Análisis Numérico”, Grupo Editorial Iberoamérica 1985. Esta obra es un clásico del Cálculo Numérico, destaca por una exposición simple y al mismo tiempo clara, con múltiples ejemplos y una descripción de los algoritmos bien diseñada. [Bo] Borse G. ”Programación en fortran 77” Anaya, 1989. En esta obra se presenta el lenguaje de programación fortran 77 con numerosas aplicaciones al análisis numérico. [Ci] Ciarlet P.G. ”Introduction à l’analyse numérique matricielle et à l’optimisation”, Masson , 1990. Con un exquisito rigor se abordan los temas básicos del Análisis Numérico Matricial y métodos de optimización, incluyendo la resolución de sistemas a través de métodos directos, iterativos y métodos tipo gradiente, así como el cálculo de autovalores y vectores propios. [Hi] Higham N. "Accuracy and Stability of Numerical Algorithms", SIAM, 1996 Esta obra, muy reciente, da una visión general sobre los últimos avances en Análisis Numérico, haciendo especial énfasis en la precisión de los algoritmos numéricos y en la propagación de errores, también resulta de interés la descripción de las aritméticas que utilizan los ordenadores más recientes como la aritmética Standard de I.E.E.E.. [Hu] Hultquist P. F. ”Numerical Methods for Engineers and Computer Scientists”, The Benjamin/Cummings Publishing Company, Inc. 1988. Esta obra, presenta una cuidada selección de temas básicos en Análisis Numérico, sin pretender ser tan exhaustiva como otras obras de carácter más general, la buena presentación de los temas elegidos la hacen de interés.
APÉNDICE A: Resumen de los comandos de UNIX En este breve resumen seguiremos el siguiente esquema. En primer lugar, aparece el comando UNIX. A continuación, entre paréntesis, su equivalente en MS-DOS (si existe). Finalmente, un comentario y un ejemplo. cd (cd) cambia el directorio activo >cd /users/p701/fortran77 more (type) visualiza el contenido de un ﬁchero >more /users/p701/fortran77/programas/prog1.f ls (dir) visualiza contenido de un directorio >ls /users/p701/fortran77 cp (copy) copia un ﬁchero en otro. >cp /users/p701/fortran77/programas/prog1.f . rm (del) borra un ﬁchero >del prog1.f man (help) suministra ayuda sobre un comando > man ls logout se termina la sesión y se sale del sistema >logout ps visualiza los números de procesos que están abiertos que corresponden al usuario alumno >ps -u alumno kill interrumpe la ejecución de un proceso de número N proceso >kill -9 N proceso mkdir (mkdir) crea un directorio >mkdir practica1 rmdir (rmdir) borra un directorio >rmdir practica1 mv (move) cambia de nombre o ubicación un archivo. >mv prog1.f practica1.f chmod cambia los permisos de lectura, escritura y ejecución de un ﬁchero. Este comando es de utilidad para salvaguardar la información de directorios y ﬁcheros de miradas ajenas. Hacer > man chmod para mirar las opciones. chown cambia el propietario de un ﬁchero. Hacer > man chown para mirar las opciones. du (tree) visualiza la cadena de directorios >du /users/p701 ﬁnd busca un archivo de nombre f ile en el directorio dir >ﬁnd dir -name f ile -print grep busca los ﬁcheros que contenga la cadena de caracteres string >grep string * APÉNDICE B: Resumen del procesador de texto vi El procesador de texto vi tiene la ventaja de estar presente en cualquier máquina que trabaje sobre UNIX y no requiere ningún entorno gráﬁco. Puede ejecutarse en dos modos. El modo comando (el que está por defecto al entrar en vi), donde se ejecutan comandos, y el modo edición, que es donde se escribe normalmente el texto.
Intercambio entre modo comando y modo edición ESC pasa de modo edición a modo comando i pasa de modo comando a modo edición A pasa a modo edición y pone el cursor al ﬁnal de la línea O inserta una nueva línea, pasa a modo edición y pone el cursor al principio de la nueva línea Manejo de Ficheros (en modo comando) :w escribe en disco el ﬁchero :wq escribe en disco el ﬁchero y sale del vi :e f ichero.name edita el ﬁchero f ichero.name :q! sale del vi sin guardar cambios. :w f ichero.name escribe el ﬁchero actual en el ﬁchero f ichero.name en disco !comando ejecuta el comando UNIX comando :set nu presenta los números de línea en pantalla Comandos para desplazarse por el texto (en modo comando) Crtl F página adelante Crtl B página atrás $ pone el cursor en el ﬁnal de la línea 0 pone el cursor en el principio de línea /string busca hacia adelante el string string ?string busca hacia atras el string string n repite la última búsqueda G va al ﬁnal del texto 3G va a la línea número 3. Comandos para borrar líneas o caracteres (en modo comando) x borra el carácter donde se encuentra el cursor r character remplaza el carácter donde se encuentra el cursor por el carácter character dd borra la línea donde se encuentra el cursor 3 dd borra 3 líneas desde donde se encuentra el cursor hacia abajo dw borra la palabra donde se encuentra el cursor Comandos para copiar y desplazar bloques (en modo comando) yy copia en el buﬀer la línea donde se encuentra el cursor 3yy copia en el buﬀer 3 líneas hacia abajo desde el cursor dd copia (y borra) al buﬀer la línea donde se encuentra el cursor 3dd copia (y borra) al buﬀer 3 líneas hacia abajo desde el cursor p copia el contenido del buﬀer en el texto.
APÉNDICE C: Algunos fallos comunes en Fortran 1. No poner EN D al ﬁnal del programa principal o de la función.
2. Escribir números como 1/2 ó 10 ∗ ∗20 en precisión entera. Solución: Escribir 1./2.ó 10. ∗ ∗20. 3. Utilizar variables enteras como ﬂotantes o al revés. Sugerencia: Aunque no sea necesario, declarar los tipos de todas las variables que se utilicen al principio del programa o función. 4. Utilizar un parámetro de una función para asignar dinámicamente memoria a un vector o matriz en el interior de la función. Solución: Poner una declaración de PARAMETER al principio de la función y con ella asignar las memorias de forma estática. 5. No poner ningún comentario en los programas. 6. Anidar excesivamente los programas. Siempre hay que buscar que el número de anidamientos sea mínimo. 7. No respetar los tipos en los pasos de parámetros de las funciones. 8. Utilizar vectores sin declararlos con la sentencia DIMENSION. 9. No pasar la dimensión de un vector como parámetro de una función. 10. Exceso de sentencias GOT O. Las sentencias GOT O pueden diﬁcultar el seguimiento del ﬂujo del programa y sólo hay que utilizarlas cuando sean indispensables. 11. A veces, los programas pueden fallar por errores de redondeo en los cálculos. Fortran da la posibilidad de cambiar el número de bits utilizados para almacenar las variables en el momento de la compilacion. Por ejemplo, si hacemos ”f77 -rn prueba.f -o prueba” donde n es 8 ó 16, aumentaremos la precisión de la aritmética para las variables reales. Análogamente, si en lugar de utilizar la directiva -rn utilizamos -dn aumentaremos la precisión de las variables declaradas DOUBLE PRECISION, y si utilizamos -in las variables enteras.
More From This UserInterpolacion de Amperaje y Voltajehelice bronceCapítulo 10 Las Reglas de Simpson ParaNautic Al ESEl Soldado de Metal Con BronceFile_2182_manual de Soldadura Indura 2007Hydra SoftwareTrabajo Externo de Calculo2Ejercicios Resueltos de Calculo Vectorial e Integrales de LineaPandeo en Columnas.ppt,Dmil14108200910Cap.1 Curvas, Superficies Funciones de Vrias Variables
Apuntes de Analisis Numerico by Francisco Peregué5,6K viewsEmbedDownloadRead on Scribd mobile: iPhone, iPad and Android.Copyright: Attribution Non-Commercial (BY-NC)Download as PDF, TXT or read online from ScribdFlag for inappropriate contentMore informationShow less

References: resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución