Source: https://es.scribd.com/document/52039575/Apuntes-de-analisis-numerico
Timestamp: 2016-10-21 22:16:25+00:00

Document:
NavegarNavegarInteresesBiography & MemoirBusiness & LeadershipFiction & LiteraturePolitics & EconomyHealth & WellnessSociety & CultureHappiness & Self-HelpMystery, Thriller & CrimeHistoryYoung AdultNavegar porLibrosAudio librosCómicsPartiturasExplorar todoSubirIniciar sesiónRegistrarseLibrosAudio librosCómicsPartiturasApuntes+de+analisis+numericoUploaded by Daniel Nuñez Hernandez11K visitaDescargaInsertarSee MoreCopyright: Attribution Non-Commercial (BY-NC)Precio de lista: $0.00Download as PDF, TXT or read online from ScribdFlag for inappropriate contentANÁ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
(TOL*(ABS(A)+10. Sea u la unidad de redondeo de dicha aritmética. Calcular el número inmediatamente inferior a él en dicha aritmética. Para evitar este comportamiento.B.A continuación. B son iguales con una tolerancia T OL (tomando el máximo de A. los criterios anteriores quedan | A |≤| A | T OL lo cual es imposible (si T OL < 1). Un resultado importante para la comparación de dos números es el siguiente: Teorema 3 Si z1 .ABS(B)) THEN IF(ABS(A-B). Un número real cualquiera z. se tiene que 2e < 2 | z | y.EQ. Por ejemplo. en una aritmética de precisión ﬁnita. Problema 9 (2 puntos) Dado un número z = e P an 2e t n=1 2n .**(10. Teorema 2 Sean ymin .A. an = 0 o an = 1. que tomaremos positivo sin pérdida de generalidad.))) THEN IGUAL=0 RETURN ELSE IGUAL=1
Ahora bien. como a0 = 1. entonces e e donde z es el número más cercano a z en la aritmética. z2 ∈ A son distintos entonces e e Demostración: Ejercicio En muchos algoritmos. También se puede utilizar un criterio más simple. con = 10−10 READ *. si B = 0.’A distinto de B segun la tolerancia TOL’ STOP ENDIF END FUNCTION IGUAL(A.B. salvo que A también sea 0. 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.**(10. | B |} + ) T OL
donde a0 = 1 y. se puede añadir al criterio un valor > 0 de la siguiente forma: | A − B |≤ (max {| A |.))) THEN IGUAL=0 RETURN ELSE IGUAL=1 RETURN ENDIF ELSE IF(ABS(A-B). B). en general. | z2 |} u e e e e
.(TOL*(ABS(B)+10.TOL). Además. 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.TOL) IF(ABS(A). se puede expresar como ∞ X an z = 2e 2n n=1
Por el problema anterior. como | A − B |≤| A | T OL pero en este caso le estamos dando una signiﬁcación especial a A con respecto a B. y por tanto | z − z |≤ e 2e 2−t 2
Programa 4 Programa en Fortran 77 que determina si dos variables A. | B |} T OL Este criterio es simétrico en el sentido de que trata de igual modo los números A y B.B. el test de parada incluye el hecho de que dos variables estén próximas entre sí.TOL IF(IGUAL(A.0) THEN PRINT *. 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 |. por tanto | z − z |≤| z | 2−t =| z | u e con lo que queda demostrado el teorema. Si un número real z veriﬁca que ymin <| z |< ymax .LE. el número que está a la derecha de la desigualdad también pertenece a la aritmética de precisión ﬁnita. ymax los valores positivos menor e e y mayor de una aritmética de precisión ﬁnita.GT. 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 |.LE.’A=B segun la tolerancia TOL’ STOP ELSE PRINT *. 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.
va a producir un pequeño redondeo. Por ejemplo. para minimizar el efecto de los redondeos en las operaciones. Nosotros no vamos a entrar en este curso en cómo se pueden deﬁnir algorítmicamente estas operaciones. el resultado sí es exactamente 1. en los algoritmos. se deberá evitar. Son los que se producen al ”redondear” un número real para poder expresarlo en una aritmética de precisión ﬁnita. que es donde más error hay. en el programa Fortran anterior. existen 4 operaciones básicas. a menudo. se cancelan las partes signiﬁcativas. para ser más precisos numéricamente.1 el ordenador. parece razonable pensar que. pero. se ha utilizado este fenómeno de cancelación para poner de maniﬁesto la diferencia entre trabajar con bases distintas. números tan naturales para nosotros como 0.2**7 B=B+A CONTINUE DO 2 K=1. Como vimos en la sección anterior.
C=0. A=2**(-7. 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 . al representar 0. el resultado sea exactamente 1. Como conclusión de este apartado.01 evitando los errores de cancelación. los errores de redondeo se van acumulando en la parte menos signiﬁcativa del número (los dígitos de menos valor). cuando sumamos 100 veces el número 0.(1-B)*(10**10) PRINT *. este programa permite identiﬁcar la base de la aritmética con la que trabaja el ordenador. a continuación. Solamente queremos mencionar que. quedando la aportación de los dígitos de menos valor. en la conocida fórmula del cálculo de raíces de un polinomio de grado 2. 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. este error está controlado por la denominada unidad de redondeo. realizar la operación en una aritmética de mayor precisión y. que corresponde a los dígitos de mayor valor.) B=0
. las aritméticas estándares de ordenador trabajan en base 2. Por ello. ﬁnalmente.01 D=0 DO 1 K=1. Estos errores se producen al restar números de aproximadamente la misma magnitud.01.100 D=D+C CONTINUE PRINT *. la resta de variables que tengan una magnitud cercana. dejando relativamente intacta la parte más signiﬁcativa del número. Por ejemplo. al realizar operaciones sobre una variable. de tal forma que. al restar dos números de magnitud parecida. Esto quiere decir que. Sin embargo. podemos extraer que. u = 2−t .RETURN ENDIF ENDIF END 1 Asociado a cualquier aritmética de precisión ﬁnita de números reales. destacaremos 3 tipos: Errores de redondeo. Sin embargo. Por ejemplo. la resta. muchas veces se intenta evitar la posibilidad de restar 2 números que pudieran ser de magnitud parecida. 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. 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. si sumamos 128 = 27 veces el número 2−7 . Este tipo de errores se produce al realizar un cambio de base para representar un número real. el error de redondeo tiene la e expresión: | z − z |≤| z | u e Errores de cambio de base.1 no pueden representarse de forma exacta en una aritmética en base 2. Como vimos en la sección anterior. la multiplicación y la división de números reales dentro de la aritmética. Problema 10 (1 punto) Calcular las raíces del polinomio P (x) = x2 − 2x + 0. Hay que tener en cuenta que. y este pequeño error de redondeo se puede ir propagando hasta producir errores apreciables. al tomar un número real z y aproximarlo en la aritmética por el valor z ∈ A más próximo. no es así. los humanos pensamos y razonamos en términos de números en base 10. En los algoritmos. que es como solemos hacerlo. el resultado se redondea para pasarlo a la precisión inicial. Por ejemplo. 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.(1-D)*(10**10) END
Fuentes de errores numéricos Dentro de las posibles fuentes de errores numéricos. Por lo tanto. en la medida de lo posible. Errores por Cancelación. que son la suma.
obtenemos. 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. 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. 2 si f ( a+b ) · f (b) < 0 entonces hacemos a = a+b y volvemos 2 2 a subdividir el nuevo intervalo [a.*X0) IF(IGUAL(X0.X1.
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]. Si. f (x).
.0) THEN PRINT *. Si f ( a+b ) = 0 ya 2 2 hemos encontrado la raíz x = a+b . b] van aproximando la raíz. f (a)) y (b. 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.
Método de la bisección Se considera un intervalo [a.Nmax X1=X0-(X0*X0-A)/(2. es decir f (x) ≈ f (x0 ) + f 0 (x0 )(x − x0 ) que corresponde a un polinomio de grado 1. f (a) · f ( a+b ) < 0. En caso contrario. de forma general.TOL).’No máximo de iterac. y un número máximo de iteraciones N max . por el contrario. b] donde la función f (x) cambia de signo. en el razonamiento anterior.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.Nmax IF(A.’LA RAIZ DE A ES’. El método consiste en ir dividiendo el intervalo [a. Es decir.
Método de Newton-Raphson Éste es. es decir f (a)·f (b) < 0. 2]. b] por la mitad de la siguiente forma: Se toma el punto medio a+b . es decir: f (x0 ) x1 = x0 − 0 f (x0 ) y por tanto. f (b)) con el eje x. Las sucesivas subdivisiones del intervalo [a.TOL. y a continuación se calcula x1 como el cero de este polinomio. los valores de x para los cuales f (x) = 0. a partir de x0 una secuencia xn de valores que van aproximando la raíz. Dada una aproximación inicial de la raíz x0 .’El numero A no es positivo’ STOP ENDIF X0=(1+A)/2.0) THEN PRINT *.LE. una aproximación mejor x1 de la raíz. entonces f (x) = x2 − A = 0.X0 STOP ELSE X0=X1 ENDIF CONTINUE PRINT *. 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. b]. READ *. esto es. 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. sin duda.EQ. uno de los métodos más importantes y útiles para el cálculo de raíces.A. 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. DO 1 K=1. a partir de x0 . se busca. entonces hacemos b = a+b y volvemos 2 2 a empezar. Problema 15 (3 puntos) Escribir el código en Fortran 77 para implementar el método de la Regula-falsi. se considera el punto de intersección de la recta que pasa por los puntos (a. de teniendo en cuenta que si x = A.
0). 0). x2 = (1. el número máximo de iteraciones N max. 0). 0). 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. La función a la que se le calculan los ceros se deﬁne en el propio cuerpo del programa. 0). x0 y x1 . IMPLICIT COMPLEX (C) CX=(-1. Utilizar el método para calcular los posibles ceros de las siguientes funciones: 1. 0). 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. 0). se sustituye el valor f 0 (xn ) en el algoritmo. Para el ejemplo 2 tomar como datos iniciales x0 = (3. 0). 0). Para el ejemplo 1 tomar como datos iniciales x0 = (3. x1 = (2. 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. Para el ejemplo 4 tomar como datos iniciales x0 = (3. En el caso en que no conozcamos analíticamente el valor de la primera y segunda derivada de f (x). Programa 7 Programa en Fortran 77 donde se muestra un ejemplo de manejo de números complejos. Para el ejemplo 3 tomar como datos iniciales x0 = (3. f (x) = (x2 + 1)x 3. f (x) = x − 2 5. x1 = (2. son necesarias dos aproximaciones iniciales. (xn−2 . En este caso. 0). f (x) = ex − 1 4. T OL. Dicha raíz será la aproximación xn de la raíz de f (x) en la etapa n. f (xn−1 )).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.0001 y N max = 100. x0 y x1. f (xn−2 )) y (xn−1 . x2 = (1. podemos utilizar las siguientes aproximaciones: f 00 (xn−1 ) ≈ 2
Nota: Utilizar como tolerancia T OL = 0. 0). En el caso en que f 00 (xn−1 ) = 0. x1 = (2.01. x1 = (0. N max. Crear un programa en Fortran 77 que tenga como datos de entrada: las tres primeras aproximaciones de la raíz. 0). como por ejemplo polinomios. 0) x2 = (1.0) CY=CF(CX) PRINT *. 0). x1 = 1. x1 y x2 . para determinar la igualdad entre dos números.
. x0 . Para obtener una aproximación xn mejor de la raíz calculamos los ceros del polinomio de segundo grado anterior. 0). Para el ejemplo 5 tomar como datos iniciales x0 = (3. y calcular posteriormente las derivadas de dicha parábola. Es una generalización del método de Newton-Raphson. y la tolerancia. el número máximo de iteraciones. 0).de una forma analítica. f (xn−3 )) . para determinar la igualdad de dos números. 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). en el sentido de que. nos quedamos con los términos hasta el orden 2. x2 = (1. la elección de las fórmulas anteriores equivale a aproximar f (x) por la parábola que pasa por los puntos (xn−3 . 0) y x0 = (1. 0). f (x) = x2 + 1 2. en lugar de quedarnos con la parte lineal del desarrollo de Taylor de la función. x1 = (2. x2 = (1. x1 = (2. 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. Los datos de entrada son las aproximaciones iniciales. y la tolerancia T OL.CY END FUNCTION CF(CX) IMPLICIT COMPLEX(C) CF=SQRT(CX) END
Práctica 2 (Método de Müller. x2 = (0. f (x) = 1
Como veremos posteriormente. 4 horas) Implementar el método de Müller. nos quedamos con aquélla que sea más cercana a xn−1 . calculamos xn por el método de Newton-Raphson.
Cálculo de las raíces de un polinomio Los polinomios son un tipo particular de funciones que.000.. éste se puede expresar también de la forma siguiente: P (x) = a0 + x (a1 + x (a2 + x (a3 + x(. El siguiente
. Ahora bien.’P(X)= ’.. Veamos que se veriﬁca que P (x) = (x − x0 )Q(x) + b0 Efectivamente. requieren el uso de alguna de las técnicas presentadas en esta asignatura. + i)30−n = 26. Sea P (x) = an xn + an−1 xn−1 + .. debemos calcular las raíces del polinomio en i dado por P (i) =
(100.. el segundo año 110.. entonces. 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. Además.. Si llamamos i al interés nominal anual que se aplica al dinero.000) (1. etc.. ’Grado Superior al Maximo’ STOP ENDIF PRINT *. es decir el primer año 100. 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 . almacenándolos en las variables P X y P P X..000
Ahora bien. Teorema 4 (Método de Horner). el método de Müller.. dado que ak = bk − bk+1 x0 y an = bn . + a0 ..487%. Polin. obtenemos P 0 (x) = (x − x0 )Q0 (x) + Q(x) de donde sale obviamente que P 0 (x0 ) = Q(x0 ).X) PRINT *. Este ejemplo muestra como un problema ﬁnanciero sencillo nos lleva a la necesidad de calcular los ceros de un polinomio.000. 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.000 pesetas todos los años. con frecuencia.. para calcular i. + (b0 − b1 x0 ) Por último.PPX END
(100.NMAX) THEN PRINT *.A.
Algoritmo de Horner para evaluar un polinomio en un punto Dado un polinomio P (x) = an xn + an−1 xn−1 + .. una vez terminada la carrera y en el desarrollo de la actividad profesional.. como muestra el siguiente programa Fortran. + b1 0 0 Demostración Sea el polinomio Q(x) = bn xn−1 + bn−1 xn−2 + .. por su gran utilidad.. + i)30−n − 26. para su resolución. + x (an−1 + xan ))))) .. aportación que se va incrementando cada año en un 10%.. aunque también hay que indicar que existen algoritmos versátiles para el cálculo de las raíces complejas.N IF(N. necesitamos evaluar tanto el polinomio como su derivada.. Nos ocuparemos sólo de las raíces reales de los polinomios. Las entidades ﬁnancieras venden a sus clientes los planes de pensiones de la siguiente forma. 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. Este teorema permite calcular el polinomio y su derivada en un punto de forma muy sencilla. como.’P ‘(X)= ’. si queremos utilizar un método de cálculo de raíces como el de Newton-Raphson. obtenemos la igualdad anterior teniendo en cuenta que (x − x0 )Q(x) + b0 = bn xn + (bn−1 − bn x0 )xn + .GT..1) (1. visto anteriormente. por ejemplo: si usted aporta durante 30 años 100. Programa 8 El siguiente programa en Fortran 77 calcula la evaluación de un polinomio y su derivada en un punto X.’Escribir Grado del Polinomio’ READ *. le aseguramos que al ﬁnal del trigésimo año tendrá a su disposición la cantidad de 26.. PARAMETER(NMAX=1000) DIMENSION A(0:NMAX) COMMON/POL/PX..N READ *. por ejemplo..X CALL HORNER(N. aparecen problemas que.. requieren un análisis algo más detallado.. El siguiente ejemplo es una buena prueba de ello.000.PX PRINT *.000
El cálculo de las raíces de este polinomio nos lleva a i = 4.000 de pesetas. A menudo.000. + a0 .1) (1.000..... Ejemplo 2 Actualmente están muy de moda los planes de pensiones.’ DO 1 K=0..’EscriBir Coef.PPX PRINT *.A(K) PRINT *..000) (1.’Escribir valor de X’ READ *. Mi experiencia como docente en esta disciplina es que.
P 00 (x). La declaración DIM EN SION A(0 : N M AX) deﬁne un vector A. debe incluir en su inicio la misma sentencia COM M ON... + a0 con an 6= 0... se aplica el teorema anterior cambiando x por −x.... Si cambiamos x por −x..1
SUBROUTINE HORNER(N.
Problema 21 (1 punto) Calcular el número máximo de raíces positivas y negativas del polinomio x5 −35x3 +30x2 + 124x − 120... entonces el número de raíces positivas es igual al número de cambios de signo en los coeﬁcientes an . Demostración Teorema de Rolle.+ak x x (n − k)! (n − k − 1)! 1
1 − |x| ≥ k=0.. etc... 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 . + a0 . En este caso... para intercalar los ceros de una derivada con los ceros de la siguiente..n−1 |ak |) = >0 |x| − 1 = |an | |x|n − max |ak |
Demostración Es inmediato. Para la estimación del número de raíces reales negativas.n−1 | ak | maxk=0. Para que una subrutina pueda hacer uso de esas variables. P P X deﬁne la zona de memoria denomina P OL donde se encuentran las variables globales P X.. |P (x)| ≥ |an xn | −
k=0. hay 3 cambios de signo y hay una o tres raíces negativas. ≤ x0 1 2 m−1 ≤ xm ≤ Pmax Volviendo a aplicar este razonamiento sucesivamente sobre P 0 (x).-1 PX=PX*X+A(K) PPX=PPX*X+PX CONTINUE PX=PX*X+A(0) END
Demostración [Is-Ke] Pg. −3.
.n−1 |ak | . los signos de los coeﬁcientes son + − +−. es decir −Pmax ≤ x1 ≤ x0 ≤ x2 ≤ x0 ≤ .. entonces las m raíces distintas |an | x1 < x2 < .n−1 | ak | −1 − .. −1.n−1 maxk=0. evaluar el polinomio y su derivada en el punto x = 2.PPX PX=A(N) PPX=A(N) DO 1 K=N-1....1..n−1 1 − |x| n |x| n = ≥ |an | |x| − max |ak | k=0. entonces las raíces reales de P (x) están en el intervalo ∙ ¸ maxk=0... . Por tanto. los signos de los coeﬁcientes son: + + −−. Teorema 8 La derivada k − esima P k) (x) del polinomio ´ P (x) = an xn + an−1 xn−1 + . Ejemplo 3 Sea P (x) = 3x4 + 10x3 − 10x − 3...X) DIMENSION A(0:*) COMMON/POL/PX. utilizando el algoritmo de Horner.... derivando sucesivamente el polinomio P (x).... de tamaño N M AX + 1. La declaración COM M ON/P OL/P X. 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 + . < xm de P (x) están intercaladas con las 0 raíces x0 < x0 < . − 1 ...1 + | an | | an | Demostración Veamos que si |x| > 1 + entonces |P (x)| > 0.. y numerados desde 0 hasta N M AX.n−1 |ak | Pmax a 1 + .. a0 (saltando los posibles coeﬁcientes nulos). o bien ese mismo número menos un número par..
Nota: La declaración P ARAM ET ER(N M AX = 1000) permite deﬁnir constantes. < x0 1 2 m−1 de P (x).. Teorema 7 Entre dos raíces de una función derivable f (x) hay una raíz de f 0 (x).. y localizarlas en un intervalo..A. Por tanto.a0 es P k) (x) = k! an n! n−k an−1 (n − 1) ! n−k−1 + +. 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 + ..... 126. Efectivamente. las raíces son x = 1... Los dos resultados anteriores permiten aislar las posibles raíces de P (x) de la forma siguiente: Si llamamos maxk=0.n−1 |x| − 1 |x|n (|an | (|x| − 1) − maxk=0.. |an |
Problema 20 (2 puntos) Dado el polinomio P (x) = 2x3 + 3x2 + 4x + 5. P P X. hay un único cambio de signo y hay una raíz positiva.... de reales en precisión simple..
25 2. 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.166]. 2. pero que utilizan técnicas más complejas.cuyas raíces son x = −1. buscamos las raíces de P 00 (x) en esos intervalos. 3.5 -1. 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. por tanto los ceros de P 00 (x) estarían en los intervalos [−2.5 -1. 2. 30! = 2. utilizando cualquier método numérico de los vistos anterioremente. habremos aislado completamente a las raíces de P (x). 574. que tiene por raices x = 1. es decir x = 0.. Se calcula la raíz x1 del Polinomio P n−1) (x) (que es un polinomio de grado 1) 3. <
Al ﬁnal del procedimiento. puesto que su utilización requiere el cálculo de factoriales.25 0 1. 1. por ejemplo.5 x -2. 6 × 1032 . 0. el método de la Regula-falsi. Para k = n − 2. obteniendo −0. 05 × 10−2 . Para este polinomio. y cuya gráﬁca es
20 50 0 -2. las raíces están en el intervalo [−8.. 1 −Pmax <
-2.858. 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.. 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). Puesto que hay cambio de signo de P 00 (x) en cada uno de estos intervalos.5 x
La derivada de este polinomio es P 0 (x) = 4x3 − 3x2 − 14x + 1.25
50 25 0 0 -25 -50 -75 -100 1.25] y [0.858 para el inter12
. Existen métodos mejores para el cálculo de raíces de polinomios. −1.25 0 1.5 x
La derivada tercera de este polinomio es P 000 (x) = 24x − 6.25 2.166. Ejemplo 4 Consideremos el polinomio P (x) = x4 − x3 − 7x2 +x+6. cuya raíz es x = 0.25.. . 7. Por tanto.5 x
< .25 -25 25 0 0 1.25 2. El método presente en el siguiente programa.cuyas raíces son x = −0.25 2.5 -1. 8]. Se parte del intervalo [−Pmax . funciona razonablemente bien para grados de polinomios pequeños. Por ejemplo.1..5 -1.25. 358 y cuya gráﬁca es
25 0 -2. −2. Este procedimiento se puede utilizar para grados relativamente pequeños (n < 30). tenemos que Pmax = 8.25. que se dispara rápidamente. Pmax ]
858].0) THEN R(M)=PI(K) M=M+1 ENDIF 8 CONTINUE ICEROPOL=M END
FUNCTION ICEROPOL(A.N 2 F(K)=F(K-1)*K *** Calculo intervalo inicial PMAX=ABS(A(0)) DO 3 K=1. −0. X2]. [−0.L-1) 6 CONTINUE 4 CONTINUE *** Pasamos las raices al vector R() M=0 DO 8 K=1.X1.TOL PRINT *. ’Grado Superior al Maximo’ STOP ENDIF PRINT *. X1. las posibles raíces de P 0 (x) estarán en los intervalos [−4. tomando como valor inicial el punto medio del intervalo [X1. 05×10−2 . [−1.5.253.N IF(R(K-1).R.PI(L1). los posibles ceros de P (x) estarán en los intervalos [−8.TOL. 8]. 05×10−2 ].GT. para NewtonRaphson’ READ *. Dicha subrutina tiene como parámetros un vector A().Nmaxx CALL HORNER(N.N-K PI(L)=RP(N-K. DO 2 K=1.858.-1 DO 5 L=0.’ DO 1 K=0.R(0:*) COMMON/POL/PX.EQ.) THEN
. 4.A.
PARAMETER(NMAX=30) DIMENSION A(0:NMAX). 7.X2. un vector R(). También se deﬁne la función auxiliar RP (N. ’Escribir Tolerancia’ READ *.Nmaxx) PRINT *.M-1 9 PRINT *.’Escribir Coef. que devuelve la raíz del polinomio que se obtiene aplicando el método de Newton-Raphson.166. T OL. para el proceso de Newton-Raphson.25] y 1. −1.25. 574]. R(0:*).Nmaxx) PARAMETER(NMAX=30)
FUNCTION RP(N.R. ﬁnalmente. T OL. y el número máximo de iteraciones N max xx. Buscamos.AP.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. la tolerancia T OL. Por tanto. PI(0)=-PMAX PI(1)=-(A(N-1)*F(N-1))/(A(N)*F(N)) DO 10 K=2. 574. A. 1.TOL. [7.LT.R(0:NMAX-1) COMMON/POL/PX.N 7 PI(K)=PMAX DO 4 K=N-2.NMAX) THEN PRINT *. el grado del polinomio N.’Escribir Grado del Polinomio’ READ *. R.PI(L).A.EQ. 3. iter.EQ. Polin. las raíces de P (x) en cada un de esos intervalos y obtenemos x = −2.N. con la que consideramos que dos números son iguales.N. donde se guardan las raíces del polinomio una vez calculadas.N 1 READ *.358.PPX PRINT *. Problema 22 (2 puntos) Aislar en intervalos las raíces del polinomio P (x) = 20x3 − 45x2 + 30x − 1.PPX **** Calculo de los factoriales F(0)=1.L) DIMENSION A(0:*). donde están los coeﬁcientes del polinomio.0.N-1 IF(PMAX. 7.574.R.’El Pol.M. 2. −1. 05 × 10−2 y 2. PI(0:NMAX+1) COMMON/POL/PX. AP(0:NMAX).R(K) END
DIMENSION A(0:*). X2. que devuelve las raíces reales de un polinomio. 2. Max.166]. obteniendo x = −1.Nmaxx.’ raices’ DO 9 K=0. N maxx.TOL.
Programa 9 Programa en Fortran 77 donde se implementa la función ICEROP OL(A. Buscamos ahora las raíces de P 0 (x) es esos intervalos.0. 0.5]. N.N IF(N. L). 253] y [2. 1. Por tanto. DO 1 K=1.Nmaxx.R.EQ. 358 para el intervalo [0.) THEN IF(PX.Nmaxx M=ICEROPOL(A.PPX R(L)=1. ’Escribir No. tiene’.valo [−2.N 10 PI(2)=PMAX *** Calculo de los coeﬁcientes del *** polinomio derivada DO 7 K=2. N maxx).358] y [1. 253. IF (X1.ABS(A(K)) THEN PMAX=ABS(A(K) ENDIF 3 CONTINUE PMAX=PMAX/ABS(A(N))+1. R.0.RP) IF (PPX.A(K) PRINT *.TOL. F(0:NMAX).X2) THEN RP=X1 RETURN ENDIF RP=(X1+X2)/2.
Es decir. 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. 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.5
0. es el único polinomio de grado menor o igual que N tal que PN (xi ) = f (xi ) ∀i = 0.TOL). lo cual es imposible. x1 = −1 y x2 = 1.0) THEN RP=RP1 R(L)=0. Por tanto Q(x) ≡ 0 y P (x) = PN (x).N .. N
-1..N . .
INTERPOLACIÓN DE FUNCIONES I El problema general de la interpolación de funciones consiste en. posee N + 1 raíces...RP. vamos a interpolarla en los puntos x0 = 0.EQ. N.1
R(L)=0. .
. 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. sabemos que f (xi ) = fi .. 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. salvo que Q(x) sea identicamente igual a cero. π . Para calcular P2 (x). aproximar el valor de la función fuera de ese conjunto ﬁnito de puntos.. El polinomio interpolador de Lagrange PN (x) de f (x) en los puntos {xi }i=0.. 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. .. 2 2
estos polinomios base tienen la propiedad fundamental siguiente ½ 1 si i = j P i (xj ) = 0 si i 6= j Por tanto..5 -1 -0. RETURN ELSE RP=RP1 ENDIF ENDIF CONTINUE END
Ejemplo 5 Consideremos una función f (x) = ex . por tanto. Entonces.
Interpolación por polinomios de Lagrange Sea una función f (x) que conocemos en un conjunto ﬁnito de valores {xi }i=0. RETURN ELSE RETURN ENDIF ELSE RP1=RP-PX/PPX IF(IGUAL(RP1. a partir del conocimiento del valor de una función (y eventualmente de sus derivadas) en un conjunto ﬁnito de puntos..5
PN (x) se puede expresar en término de los denominados polinomios base de Lagrange P i (x). el polinomio interpolador de Lagrange en estos puntos. π y 3π ..
Además. 1]. Demostración Si x = xi . Un método más directo para el cálculo de PN (x) es el denominado método de diferencias de
. Teorema 11 Sea N ≥ 0. b] se comete. debemos cambiar todos los polinomios base de Lagrange. un error de interpolación. Por tanto. 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. los valores óptimos de interpolación xi dados por la fórmula anterior son las raíces de los denominados polinomios de Chebychev.Error de interpolación de Lagrange y polinomios de Chebychev Evidentemente. b] Se consideran los puntos xi dados por µ µ ¶¶ b−a 2i + 1 1 + cos π i = 0.b] f N +1) (ξ) (N + 1)!2N µ b−a 2 ¶N+1
para cualquier otra elección posible de valores de interpolación xj . N xi = a + 2 2N + 2
Ejemplo 6 Se considera [a. si queremos añadir un nuevo punto de interpolación. 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. 1] se encuentra en [Ki-Ch] Pg. 2π] interpolando en los puntos 0. entonces f (x) − PN (x) = f N +1) (ξ) N Π (x − xi ) (N + 1)! i=0
x∈[a. π y 2 3π 2 . 629 41 . π . 1] al interpolar la función cos(x) en los puntos descritos en el ejemplo anterior. 703 7 × 10−2
Problema 26 (2 puntos) Calcular el error máximo de interpolación en el intervalo [0. La cuestión que vamos a abordar en este apartado es. y que nosotros podamos elegir los valores de interpolación xi . Por tanto. 146 45 1. b]. y un intervalo [a. al aproximar f (x) por el polinomio interpolador PN (x) en un intervalo [a. el error de interpolación máximo viene determinado por: | f (x) − PN (x) |≤ maxx∈[a.b]
Demostración La demostración para el intervalo [−1.. despejando y sustituyendo λ por su valor. obtenemos el resultado del Teorema. y PN (x) su polinomio interpolador de Lagrange en los puntos {xi }i=0. La demostración para un intervalo cualquiera [a. b]. elegiremos los puntos xi tales que ΠN (x − xi ) sea lo más pequeño posible en i=0 [a. Para ello. b]. 1] y N = 5 (es decir 6 puntos de interpolación). φ00 (t) tiene al menos n − 1 ceros y así sucesivamente hasta llegar a φN+1 (t).. 1] en [a.N ⊂ [a. 982 96 .. en el caso en que queramos interpolar una función en un intervalo [a. b] se obtiene fácilmente transformando el intervalo [−1. cómo elegirlos de tal forma que el error de interpolación sea mínimo.. 370 59 . y acotarlo superiormente. Si llamamos ξ a dicho cero. 853 55 . b] = [−1. 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. En el caso de que [a. Análogamente. Los puntos de interpolación dados por el teorema anterior son: x0 x1 x2 x3 x4 x5 = = = = = = . utilizando este resultado. su función derivada φ0 (t) tiene al menos n ceros repartidos entre los ceros de φ(t). el error de interpolación es cero y por tanto la fórmula anterior es válida.b]
x∈[a. b]. 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. que tiene al menos 1 cero. e
donde ξ es un valor intermedio perteneciente a [a. b] = [0. b]. . 292-294. en general. b] y x ∈ [a.. TN (x). obtenemos φN+1 (ξ) = f N+1) (ξ) − λ(N + 1)! de donde. Teorema 10 Sea f (x) una función.. que viene determinado por el siguiente teorema.
los coeﬁcientes f [x0 .Newton. 2] = 2 e3 − 2e2 + e1 f [1. 1. .. . . veriﬁcan las siguientes propiedades: f [xi ] = f (xi ) f [xi+1 ] − f [xi ] f [xi .. .. xk de la siguiente forma: P0 (x) = a0 P1 (x) = P0 (x) + a1 (x − x0 ) P2 (x) = P1 (x) + a2 (x − x0 )(x − x1 ) . . xi+k−1 ] f [xi . xi+1 ] = f [xi+1 . P0 (x) = 1 P1 (x) = 1 + a1 x Como P1 (1) debe ser igual a e. xk ] Ejemplo 7 Vamos a interpolar la función f (x) = ex en los puntos x0 = 0... por tanto: f [xi . ... . xi+k ] = f [xi+k ... . xi+k ] indica...... observamos que f [xi . . xi+k ] = f [xi+k . xi+k ] − f [xi . . .... donde bj = f [xi+k . x2 = 2.. xk ].. .. obtenemos el polinomio interpolador de la siguiente forma: f [0... xi ] = bk De nuevo. los coeﬁcientes que acompañan a la potencia xk−1 en ambos polinomios coinciden y... .. despejando obtenemos a1 = e − 1 Por último P2 (x) = P1 (x) + a2 x(x − 1)
Demostración En primer lugar... xi+k y..... 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... xi+k−j ] Por la unicidad del polinomio interpolador obtenemos que Pk (x) = Qk (x) y. xi ] Consideremos ahora el polinomio interpolador Qk (x) que interpola en los puntos xi+k ..... 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..... para cada Pk (x)... 2. xi+k ] no depende del orden en que tomemos los puntos xi . que se denominan diferencias divididas de Newton. xk ] veriﬁcan f [xi+1 .. 2] = e2 − e1 f [2...... PN (x) = PN −1 (x) + aN (x − x0 )(x − x1 )... xi+k−1 ] f [xi .. . x1 = 1 y x2 = 2. 1] = e1 − 1 f [1. f [xi ... xi .. x1 = 1. el coeﬁciente que acompaña a la potencia xk en el polinomio interpolador Pk (x) para los puntos xi . xi+k ] − f [xi ... . .. Qk (x) se puede escribir como Qk (x) = b0 + b1 (x − xi+k ) + b2 (x − xk+i )(x − xk+i−1 ) + ... cambiando el orden de los puntos... por tanto ak = f [xi ... . xk ]...(x − xN−1 ) A los coeﬁcientes ak los denotamos por ak = f [x0 . .. teniendo en cuenta que ak−1 bk−1 = f [xi . El método consiste en ir calculando progresivamente los polinomios Pk (x) que interpolan la función en los puntos x0 .. xi+k ] = xi+k − xi
... xi+k ..... . f [xi+1 ... . 3] = 2 e3 − 3e2 + 3e1 − 1 f [0.. 2. xi+k ]
Como veremos en el teorema siguiente.. xi+1 ] = xi+1 − xi . 3] = e3 − e2 e2 − 2e + 1 f [0. despejando obtenemos e2 − P1 (2) a2 = 2 Por tanto. Como el polinomio interpolador es único.. xi+k−1 ] = f [xi+k . por la unicidad del polinomio interpolador.. es decir. 1. . x3 = 3. 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 .. por tanto: ak−1 − ak
Como P2 (2) debe ser igual a e2 . . xi+k ] = xi+k − xi Teorema 12 Si denotamos por ak = f [x0 .
Expresar el polinomio tomando en primer lugar x0 = 0. π y 3π utilizando las diferencias divididas 2 2 de Newton.EQ.0. 3 y 4 utilizando las diferencias divididas de Newton.N N=N-1 PRINT *.5 1 1. 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 ).’Test de Comprobacion’ DO 4 K=0. x2 = 1.N READ *.N) Parameter(Nmax=1000) DIMENSION A(0:*). Problema 28 (2 puntos) Calcular el polinomio interpolador de Lagrange P3 (x) de la función f (x) = sen(x) en los puntos 0.B(0:Nmax) DO 1 K=0.X(0:1000) PRINT *.1 -0.N PRINT *.EQ. devuelve el valor de la evaluación del polinomio de Lagrange en el punto X0.N-K IF (X(K+L).F.F(0:1000).F(K). 1.’ READ *. x2 = 3 y x3 = 4 y. x1 .N READ *.’Coef.X0. en segundo lugar.N DO 3 L=0.5 x 3
3 Problema 27 (2 puntos) Interpolar la función f (x) = 10 x2 +1 en los puntos x0 = −2.’Introducir No. 3] :
PARAMETER(Nmax=1000) DIMENSION A(0:1000).1 0.X(K).F.N PRINT *. x1 = −1. π . x0 = 4.
FUNCTION EVDIFNEWTON(A.EVDIFNEWTON(A. devuelve el vector A(0 : N ) de coeﬁcientes de diferencias divididas que deﬁnen el polinomio de Lagrange (A(K) = f [x0 . f (xn−1 )) .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.X.X.N) DIMENSION A(0:*).’Introducir Valores de F()’ DO 2 K=0..X(0:*) EVDIFNEWTON=A(N) DO 1 K=N-1.N) que a partir de los coeﬁcientes dados por el vector A(0 : N ) y el conjunto de puntos de interpolación. (xn−2 . x2 = 1 y x3 = 0. Problema 30 (3 puntos) Dada una función f (x) y una secuencia de valores xn .1) THEN PRINT *. x1 = 1.5 2 2. Polinomio’ DO 3 K=0. xK ] ). f (xn−2 )) y (xn−3 . f (xn−3 )).N B(K)=F(K) A(0)=F(0) DO 2 K=1..X.F(0:*). Programa 10 Programa en Fortran 77 donde se deﬁnen las funciones IDIFNEWTON.05 0 0 -0.N).X.’Puntos de Interpolacion repetidos’ STOP ENDIF PRINT *. Ptos Interp.-1 1 EVDIFNEWTON=EVDIFNEWTON*(X0X(K))+A(K) END
. y la función EVDIFNEWTON(A.En la siguiente gráﬁca se muestra la diferencia ex − P3 (x) en el intervalo [0.05 -0. . x1 = 3.15 -0.2 -0.X(K).X(0:*).F(K) IF(IDIFNEWTON(A. aproximar f (x) por la parábola que pasa por los puntos (xn−1 .X.N) END FUNCTION IDIFNEWTON(A.’Introducir Ptos Interpol.A(K) PRINT *.X0. X(K) PRINT *. x3 = 2 utilizando las diferencias de Newton y evaluar el polinomio en x = 0 utilizando el algoritmo de Horner.’ DO 1 K=0.25 0. 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.
como son: la raíz √ cuadrada x. las funciones elementales que todos usamos. 1. 1] en una aritmética de 32 bits. del mismo orden que la unidad de redondeo u en una aritmética de 32 bits. cuya unidad de redondeo u es del orden de 10−8 . π ] y a partir de ellas deﬁnir las funciones para 4 cualquier valor x (en radianes).π] (2π − x) si x > π sen[0. Podemos deﬁnir en4 tonces las siguientes funciones: ½ cos[0. la función ln(x). 1].π] (x) = ½ sen[0. la función ex . 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.). es decir 7 puntos de interpolación. que evalua el polinomio interpolador en un punto. 1].π] (x) = cos[0. Aproximación de la exponencial ex Un número real x siempre se puede expresar como x = m + x0 . Efectivamente. lo que simpliﬁca el cálculo numérico. π ]. 2 horas) Crear una función en Fortran 77 que devuelva el valor de ex con x ∈ [0.. 0.2π] (x) = sen[0. evaluar y mostrar P6 (x) y ex .
. π ] (x) a las funciones trigonométri4 4 cas deﬁnidas sobre el intervalo [0. utilizando algunas relaciones trigonométricas es suﬁciente deﬁnir las funciones cos(x) y sen(x) en el intervalo [0.π] (x) si x ≤ π −sen[0. π ] ( 2 − x) si x > π 4 4 sen[0. π ] (x) si x ≤ 4 π π ( cos[0.. Nota: Utilizar las funciones de an. a partir de estas operaciones. . el error relativo es menor que 6. etc.. Asi.). obtenemos ya la mejor aproximación posible de ex en el intervalo [0. tomando un polinomio de grado N = 6. las funciones trigonométricas: sen(x) cos(x) y tan(x). − x2 x4 x2n + + . los desarrollos de Taylor y los algoritmos de búsqueda de ceros (como √ vimos anteriormente para x). que calcula el polinomio interpolador a partir de los puntos y valores de interpolación. resta. de em . Puesto que estas funciones son 2π periódicas. Introducir por teclado un valor x. π ]. Dado que 0 ex = em ex podemos descomponer el cálculo de ex en el cálculo. π ] (x) si x ≤ 2 sen[0. división). π ] (x) y sen[0. El error relativo es | Pn (x) − cos(x) | (x)2n+1 ≤ tan(x) cos(x) (2n + 1)! Además.
i = 0. + (−1)n 2 4! (2n)!
y el error máximo cometido por el desarrollo de Taylor en un punto x ∈ [0.. y por otro.h IDIFNEWTON(. 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).. π ] (x) si x ≤ 2 − cos[0. π ] (x) = 2 cos[0. 8 × 10−9 4 (2 ∗ 5 + 1)! Por tanto. π ] (x) = π 2 sen[0. como tan(x) es creciente en [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. donde m es un número entero y x0 ∈ [0. π ] (π − x) si x > 2
cos[0.π] (x) si x ≤ π cos[0. multiplicación. Por ejemplo. 4 ] 2 − x) si x >
sen[0. Práctica 3 (Aproximación de ex . si trabajamos con una aritmética de 32 bits. 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. 6
Comprobar que el polinomio esta bien construido. 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.6 × 10−8 y.2π] (x) =
cos[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). la función xy ..5. 2 y 3. Utilizar x = 0. π ] (x) si x ≤ π 4 4 cos[0. Las técnicas elementales para deﬁnir estas funciones consisten en utilizar la interpolación polinómica.. por un lado. en el cálculo 0 de ex para x0 ∈ [0.π] (2π − x) si x > π
El desarrollo en Serie de Taylor centrado en 0 del cos(x) es: cos(x) u Pn (x) = 1. y EVDIFNEWTON(. Utilizando como puntos de interpolación los asociados a los polinomios de Chebychev: µ µ ¶¶ 1 2i + 1 xi = 1 + cos π i = 0. denotemos por cos[0. π ] (π − x) si x > 2 sen[0. por tanto. es decir que P6 (xi ) = exi para todos los xi .. .Implementación de funciones elementales Una vez deﬁnida una aritmética en precisión ﬁnita y las 4 operaciones básicas (suma. el valor 4 máximo del error se encuentra en x = π .. es necesario deﬁnir.
en una aritmética de 32 bits tendríamos la mejor aproximación posible de la función ln(x). . Una vez obtenidos la matriz A0 y el vector b0 . no ocurre con otros métodos más rápidos. algo que. pero que requieren. por ejemplo. b = (bi ) es un vector de tamaño N que determina los términos independientes. para que el sistema sea equivalente.. y u = (ui ) es el vector solución buscado. 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. 1]. estudiaremos los métodos directos clásicos para la resolución de un sistema de ecuaciones de la forma Au = b donde A = (ai.. que para el intervalo [ 1 . por tanto. N 2 4 2N + 2 Dado que ln(1) = 0. 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. utilizando relaciones trigonométricas. N.. que la matriz sea simétrica o deﬁnida positiva. 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) ∼ ..Problema 31 (3 puntos) Aproximar la función sen(x) en el intervalo [0. Además se tiene que en el intervalo [ 1 . . 1
l=k+1 a0 k.j ) es una matriz de N xN. Para obtener A0 y b0 se calcula.l
Problema 35 (2 puntos) Calcular el número de operaciones básicas (sumas. denominado pivote. Se hace lo mismo para el vector
. Método de Gauss Este método...6931471806). Por tanto. 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. 8
ANÁLISIS NUMÉRICO MATRICIAL I En esta primera sección dedicada a la resolución de sistemas de ecuaciones lineales. el cálculo de la solución u es inmediata. tiene la gran ventaja de que se puede aplicar a todo tipo de matrices. 8
Problema 33 (3 puntos) Calcular los polinomios necesarios para interpolar las funciones trigonométricas cos(x) y sen(x) en el intervalo [0. . restas. 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. π ]. para minimizar el error relativo añadiremos como punto interpolante xN +1 = 1. π ] en una aritmética de 32 bits. es decir. utilizando las funciones ex y ln(x)?
Problema 32 (2 puntos) Demostrar que. 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. multiplicaciones y divisiones) necesarias para realizar un remonte como el presentado arriba en función de la dimensión N . 973 6 × 10−8 . el valor máximo en valor absoluto de la primera columna de A. a1 = 1 y para´ > 1 an = 0 ³P t an ó an = 1.
Para N = 10 el error máximo es 3. donde x e y son números reales. que tenga valores nulos de la diagonal hacia abajo. y se hace lo mismo con el vector b. 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 .N k = N − 1. como veremos en el futuro. A continuación. = podemos reducir el cálculo del ln(x) al rango de valores 1 2 ≤ x ≤ 1.. 1] son: 2 µ µ ¶¶ 1 1 2i + 1 xi = + 1 + cos π i = 0. en primer lugar. Problema 34 (1 punto) ¿Cómo se puede obtener la función y x . es posible calcular las funciones sen(x) y cos(x) para cualquier x (en radianes).k
a0 ul k. Utilizaremos los puntos de interpolación generados por los polinomios de Chebychev.. el número es mayor o n=1 2n igual que 1 y menor que 1. A continuación. aunque no es de los más rápidos. π ] 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. que es menor que la unidad de redondeo u y. utilizando únicamente su valor en el intervalo [0. se intercambia la primera ﬁla de A con la ﬁla donde se encuentra el pivote.
B(K) CONTINUE ELSE IF(M.B(K) M=IGAUSS(A.0) THEN PRINT *.N READ *.L) PRINT *.’SOLUCION’ DO 1 K=1. donde A es una matriz triangular inferior.N READ *. La función devuelve 0 si termina correctamente y 1 en caso contrario. N max) que resuelve un sistema.’DIMENSION DEL SISTEMA MAYOR DE LA PERMITIDA’ RETURN ENDIF DO 1 K=1. 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.b.N)=0’ ENDIF END
y el remonte da como solución u1 = 6. Parameter(Nmax=1000) DIMENSION A(Nmax.*).B(Nmax).GT .N.Nmax) PARAMETER (NmaxGAUSS=1000) DIMENSION A(Nmax. b es el vector de términos independientes. la dimensión del sistema y la dimensión máxima admitida.K)) M=K DO 3 L=K+1. La función devuelve un valor entero M que indica si se ha terminado correctamente (M = 0) o incorrectamente (M = 1.’Una columna es toda cero’ ELSE PRINT *.Nrow. u3 = 8.EQ.N.EQ.Nrow(*) DIMENSION U(NmaxGAUSS) IF (N.N PRINT *. y con ello habremos obtenido un sistema equivalente tal que la primera columna es cero de la diagonal hacia abajo.NmaxGAUSS) THEN IGAUSS=3 PRINT *.Nrow.Nrow(Nmax) PRINT *.N
.N PRINT *. Nota Importante: Las líneas de código tienen que ir todas numeradas y no pueden superar las 15 líneas. forma siguiente ⎛ ⎞ ⎛ 0 ⎝ 24 ⎠ → ⎝ 8
Programa 11 Programa en fortran 77 que implementa el método de Gauss. ’Introducir Dimension’ READ *. b.
FUNCTION IGAUSS(A.B. N. u el vector solución.1)THEN PRINT *. el vector independiente b. 2).A(K.B(*). Se deﬁne una función IGAU SS que tiene como parámetros la matriz A.Nmax).B.’Introducir vector’ DO 3 K=1. la solución se devuelve en el propio vector b. u.Nmax) IF(M.’Introducir matriz’ DO 2 K=1. y así sucesivamente hasta llegar a la mencionada matriz A0 . un vector auxiliar N row.N-1 XMax=ABS(A(Nrow(K).N DO 2 L=1. Se considera el sistema ⎛ ⎞⎛ ⎞ ⎛ ⎞ −2 −2 0 u1 0 ⎝ 6 18 12 ⎠ ⎝ u2 ⎠ = ⎝ 24 ⎠ 3 11 7 8 u3
de la misma forma. En el caso en que se ha terminado correctamente. Volvemos ahora a hacer lo mismo para convertir la segunda columna cero de la diagonal para abajo. teniendo en cuenta la siguiente relación:
Problema 38 (2 puntos) Implementar en FORTRAN la funcion IDESCEN SO(A.’A(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. u2 = −6.
Ejemplo 9 Ejemplo de descomposición según el método de Gauss.N 1 Nrow(K)=K DO 2 K=1.
Ahora bien.(2.1 ⎜ b2.2 0 0 b3. porque los errores de redondeo y de cálculo producen que esta estimación no sea exacta. Teorema 13 Sea A una matriz simétrica.IF(ABS(A(Nrow(L).LT. bn.K) DO 5 M=K.n ⎞ ⎟ ⎟ ⎟ ⎟ ⎠
Problema 40 (2 puntos) Demostrar que si A = B · B t (B triangular inferior) y |B| 6= 0. 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.i aj. El siguiente teorema da 3 posibles deﬁniciones equivalentes de una matriz deﬁnida positiva.M)C*A(Nrow(K).-1 C=0 DO 7 L=K+1.N C=C+A(Nrow(K). bn..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. . donde e es el vector e = Au − b : ErrorSistema = 1 X |ei | N |bi | + 1
Para j = i + 1. esto no suele suceder.M)=A(Nrow(L). . 0 . N 1.1.N)). un vector de términos independientes b y un vector solución u.K)) M=L ENDIF 3 CONTINUE IF(XMax.. calculado utilizando alguna técnica numérica. En el denominador se añade 1 para evitar las posibles divisiones por 0.1 ⎜ ⎝ . bn.M) 5 CONTINUE B(Nrow(L))=B(Nrow(L))-C*B(Nrow(K)) 4 CONTINUE 2 CONTINUE IF(ABS(A(Nrow(N). .N) DO 6 K=N-1.K)/A(Nrow(K).1 ⎜ B = ⎜ b3. ³ ´ P bi.
Método de Cholesky Este método sólo se puede aplicar a matrices simétricas y deﬁnidas positivas.i − i−1 b2 k=1 i.LT..i = bi.M) THEN MP=Nrow(K) Nrow(K)=Nrow(M) Nrow(M)=MP ENDIF DO 4 L=K+1.. si la solución es perfecta entonces Au − b = 0. (ii) Todos los autovalores λ de A son positivos. .L)*U(L) 7 CONTINUE U(K)=(B(Nrow(K))-C)/A(Nrow(K).i − i−1 bj.2 b3. 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.NE.1 triangular inferior.XMax) THEN XMax=ABS(A(Nrow(L).3 . Para estimar el error cometido al resolver el sistema utilizaremos la expresión siguiente. . haremos lo siguiente: dada una matriz A. N ³ ´ P 1 bj.. 0 .**(-100))) THEN IGAUSS=1 RETURN ENDIF IF(K.k k=1 Fin Para j Fin Para i
.i = ai.K) 6 CONTINUE DO I=1. mejor aproximada estará la solución del sistema.N A(Nrow(L). el algoritmo para calcular B es el siguiente Para i =r .2 . 0 b2.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. bn.3 .k bi.**(-100))) THEN IGAUSS=2 RETURN ENDIF U(N)=B(Nrow(N))/A(Nrow(N). Cuanto más pequeño sea ErrorSistema. (iii) Los determinantes de todos los menores principales de A son positivos.. entonces A es simétrica y deﬁnida positiva. 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.GT.K)).(2.N C=A(Nrow(L).
DU. ⎝ −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. para evitar tener que almacenar dos matrices.DU. 0 ⎠ ⎝ 0 . 0 ⎟ ⎜ ⎟= ⎝ 0 . simétrica.es.
Práctica 4 (Método de Cholesky.Nmax): Resuelve un sistema triangular superior.B. Nota: Normalmente. Estos archivos ejemplo sólo se pueden utilizar al compilar el programa en serdis bajo UNIX.DU. bN −1 ⎠ 0 . 0 1
Los vectores mi . 0 ⎟⎜ 0 1 . B es el término independiente y DZ es el vector donde devuelve la solución.
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. ⎝ 1 5 6 ⎠ ⎝ y ⎠ = ⎝ 12 ⎠ 4 6 17 z 36 ⎛ ⎞⎛ ⎞ ⎛ ⎞ 2 −1 0 x 1 4.DZ. el programa no reconocerá el formato de los archivos. donde DB es la matriz. Efectivamente. . mN−1 lN 0 . Hay que resolver los sistemas ejemplo y calcular el ErrorSistema.dis.N. 0 ⎜ c1 a2 .N.h. y ui se calculan utilizando el esquema: l1 = a1 1 u1 = b1 l Para i = 2.Nmax): Calcula la descomposición de Cholesky de A y la devuelve en la matriz DB.Nmax): Resuelve un sistema triangular inferior. donde DB es la matriz. Nota: El programa debe permitir introducir el sistema directamente por teclado o desde disco duro.. li. • FUNCTION ICHOLESKY(A. uN −1 ⎠ 0 .
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). . ⎝ 1 1 4 ⎠ ⎝ y ⎠ = ⎝ 12 ⎠ 4 4 17 z 36 ⎛ ⎞⎛ ⎞ ⎛ ⎞ 1 1 4 x 6 3.DZ. 0 0 l1 0 ⎜ m1 l2 0 ⎟ . Vamos a descomponer A en el producto de dos matrices triangulares de la forma siguiente: ⎛ ⎞ a1 b1 . escribiendo en la parte triangular superior de B la parte correspondiente a B t . Devuelve 0 si termina bien y −1 en caso contrario • FUNCTION ERROR_SISTEMA(A. . 6 horas) Implementar en Fortran 77 las siguientes funciones : • FUNCTION ICHOLESKY_FACTORIZACION (A. Los 3 sistemas corresponden a matrices simétricas y deﬁnidas positivas. Problema 42 (2 puntos) Calcular el número de operaciones necesarias para resolver un sistema por el método de Cholesky. Devuelve 0 si termina bien y −1 en caso contrario. basta descomponer el sistema de la siguiente forma: Bz = b Bu = z
Ambos sistemas se resuelvan rápidamente haciendo un remonte y un descenso. En este directorio hay tres ejemplos de sistemas de dimensión 10. cN −1 aN ⎞ ⎞⎛ ⎛ 1 u1 ..El interés de descomponer una matriz A por el método de Cholesky es que.Nmax): Devuelve el error cometido al resolver el sistema dado por la expresión ErrorSistema de la sección anterior.ulpgc. tanto en doble como en simple precisión. utilizando las funciones deﬁnidas en an. una para B y otra para B t . ⎟ ⎟⎜ ⎜ ⎝ 0 . • FUNCTION IDESCENSO(DB. . Devuelve 0 si termina bien y −1 en caso contrario. a continuación. 100 y 500. Devuelve 0 si termina bien y −1 en caso contrario.N. • FUNCTION IREMONTE(DB.N. Si utilizamos Linux. Resolver los siguientes sistemas ejemplo: ⎞⎛ ⎞ ⎛ ⎞ ⎛ x 6 1 1 4 1. DZ es el término independiente y DU es el vector donde devuelve la solución.B.B. es muy sencillo resolver el sistema de ecuaciones Au = b.DB. .Nmax): Resuelve un sistema por el método de Cholesky y devuelve la solución en DU.
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.N. se almacena todo en una única matriz B. N − 1
. ⎝ 1 5 6 ⎠ ⎝ y ⎠ = ⎝ 12 ⎠ 4 6 26 z 36 ⎛ ⎞⎛ ⎞ ⎛ ⎞ 1 1 4 x 6 2.
form=’UNFORMATTED’.Ndimension) CHARACTER * (*) String
***** Procedimiento que escribe en el ﬁchero String la matriz A ***** cuadrada de dimension Ndimension SUBROUTINE EscribirMatriz(String.Ndimension DO 1 K=1.A. form=’UNFORMATTED’) READ(1) A Ndimension=INT(A) PRINT *.FILE=String.Ndimension READ(1) Vector(K) 1 CONTINUE CLOSE(1) LeerVector=Ndimension END
****** Función que lee una matriz cuadrada de disco de nombre String.Vector.
Problema 45 (2 puntos) Resolver utilizando el método de Crout el siguiente sistema de ecuaciones: ⎛ ⎞⎛ ⎞ ⎛ ⎞ 2 4 0 x 6 ⎝ −1 0 4 ⎠ ⎝ y ⎠ = ⎝ 3 ⎠ 0 −1 0 z −1 Problema 46 (2 puntos) Calcular el número de operaciones necesarias para resolver un sistema tridiagonal por el método de Crout.*) OPEN(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. form=’UNFORMATTED’.Vector) CHARACTER * (*) String DIMENSION Vector(*) OPEN(1. ERR=1) GOTO 2 1 OPEN(1.FILE=String.*) OPEN(1.Ndimension WRITE(1) Vector(K) 3 CONTINUE CLOSE(1) END
STATUS=’NEW’. STATUS=’OLD’. STATUS=’OLD’.ERR=1) GOTO 2 1 OPEN(1.Ndimension WRITE(1) A(K.’DIMENSION DE LA MATRIZ ’.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.String. lo almacena ****** en la tabla Vector y devuelve la dimension del vector FUNCTION LeerVector(String.’=’. STATUS=’NEW’.Ndimension. form=’UNFORMATTED’) READ(1) B Ndimension=INT(B) PRINT *. ****** Función que lee un vector de disco de nombre String.STATUS=’OLD’.FILE=String.L) 4 CONTINUE 3 CONTINUE CLOSE(1)
.FILE=String. String.’DIMENSION DEL VECTOR ’.’=’.FILE=String.Nmax) CHARACTER * (*) String DIMENSION A(Nmax. form=’UNFORMATTED’) 2 A=Ndimension WRITE(1) A DO 3 K=1.
STATUS=’OLD’. FILE=String.
DIMENSION Vector(*) OPEN(1. lo almacena ****** en la tabla A y devuelve la dimension de la matriz FUNCTION LeerMatriz(String.Ndimension DO 2 L=1.Ndimension DO 4 L=1.
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.A. Ndimension DO 1 K=1.Ndimension READ(1) A(K. form=’UNFORMATTED’) 2 B=Ndimension WRITE(1) B DO 3 K=1.Nmax) CHARACTER * (*) String DIMENSION A(Nmax.
V) PRINT *.. obtenemos la siguiente expresión: f 0 (xi ) ≈ f (xj ) − f (xi ) + O (|xj − xi |) xj − xi
donde O (|xj − xi |) indica.Ndimension DO 4 L=1.Nmax) String=’vector. obtenemos la aproximación f 0 (1) ≈ 23 − 13 =7 2−1
Si tomamos ahora xj = 1. en el cuerpo del programa.’N=’. 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 ) + .13 − 13 = 3.Ndimension DO 1 K=1.Ndimension DO 6 L=1.N DO 2 k=1.dat’. Por lo tanto.A. utilizando el valor de f (x) en otros puntos vecinos a xi .h que se encuentra en el directorio de trabajo. Por otro lado. entonces la derivada se calcula hacia adelante.A(L.END
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. b].A(L.dat’ Ndimension=LeerVector(String. 1) (utilizar h = 0.String DIMENSION V(Nmax). diremos que el orden de aproximación es 1. String deﬁne un string de caracteres de tamaño 10.Ndimension DO 5 K=1.A.K) CONTINUE PRINT * CONTINUE END
Si tomamos un punto xj 6= xi .V(K) V(K)=2*V(K) CONTINUE CALL EscribirVector(’vector3. en este caso. b] a partir de la evaluación de f (x) en algunos puntos incluidos en el intervalo [a.
Problema 47 (2 puntos) Calcular analítica y numéricamente la matriz gradiente en el punto (1.V(K) CONTINUE Ndimension=3 DO 3 K=1. + (x − xi )N + .1) de la función: ½ 2 x + y2 − 1 f (x.Nmax) Ndimension=LeerMatriz(’matriz.K) CONTINUE PRINT * CONTINUE CALL EscribirMatriz(’matriz. mientras que si xj < xi . 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. y) = x−y
.h’ PARAMETER(Nmax=100) CHARACTER * 10. Se denomina orden de la aproximación a la potencia más pequeña que aparece en el término del error.V) PRINT *. Si xj > xi .Ndimension PRINT *.Ndimension A(L. truncamos el desarrollo de Taylor y despejamos. básicamente..N PRINT *. Consideremos la función f (x) = x3 .dat’.Ndimension PRINT *.h0 incluye el ﬁchero an.Ndimension.V. mejor será el valor aproximado de la derivada. INCLUDE ’an. la derivada se calcula hacia atrás. cuanto más próximo esté el punto xj al punto xi .1. A(Nmax.Ndimension) N=LeerVector(’vector3.. una fórmula de integración numérica es un procedimiento que permite aproximar el valor de la integral en un intervalo [a.K)=L+K PRINT *. Ejemplo 10 Veremos en este ejemplo como.. que el error cometido es una suma de potencias de |xj − xi | en la que la potencia más pequeña es 1.dat’.dat’. La declaración CHARACT ER ∗ 10. 1! N!
que está mucho más próximo al valor real.1 − 1
Nota la declaración IN CLU DE 0an. 31 1. obtenemos f 0 (1) ≈ 1.Nmax) PRINT *.
calcular el polinomio de Lagrange que interpola a f (x) en esos 3 puntos. xr . cuanto mayor es el orden de una fórmula de aproximación. f (xi − h). y) = F − hFx + 3.Problema 48 (3 puntos) Dados 3 puntos distintos xl . 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. f (xi − 2h).y) ∂x∂y . F (x + h.93 f 0 (1) ≈ = 3. y − l) = F + hFx − lFy + 1 (h2 Fxx − 2hlFxy + 2 ¡ ¢3 l2 Fyy ) + O( h2 + l2 2 ) 8. F (x. 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. Utilizaremos la notación Fi. calcular el polinomio de Lagrange que interpola a f (x) en esos 3 puntos. el orden de aproximación es 2.y) . 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. y). F (x + h. comprobamos que. y) = F + hFx + 2. calcular la derivada segunda de ese polinomio en xi . Fy =
∂F (x. F (x. y demostrar que. f (1) = 0. y + l) = F + lFy + 4.
1.1 1.y) ∂y2
Si tomamos ahora xj = 0.y) . se utilizan los desarrollos de Taylor siguientes en 2 variables. F (x − h. En general. y − l) = F − lFy +
6. f (xi + h). y − l) = F − hFx − lFy + 1 (h2 Fxx + 2hlFxy + 2 ¢3 ¡ l2 Fyy ) + O( h2 + l2 2 ) 7. f (4) = 9
aproxima la derivada de f 0 (xi ) con un orden de aproximación de 2. Fy (x. que es el vector de derivadas parciales. Problema 53 (3 puntos) Dados 3 puntos xl < xi < xr . Nótese que. la expresión anterior nos da f 0 (1) ≈ 23 − 03 =4 2
Diferenciación numérica en dimensiones superiores Estudiaremos. Deducir como queda la fórmula anterior para aproximar la derivada segunda. y comprobar que da la misma fórmula que la presentada en el problema anterior. xi . 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.y) Utilizaremos la siguiente nomenclatura: Fx = ∂F∂x . la precisión es mayor que con la fórmula anterior. y xr = xi + h. 01 0. Problema 51 (3 puntos) Dados 3 puntos. la aproximación de las derivadas de una función de varias variables. ∂y ∂ 2 F (x. utilizando f (xi ). xr . Para simpliﬁcar la exposición. (x. F (x − h. Ejemplo 11 Veremos en este ejemplo como. ∂x2
∂ 2 F (x. y xl = xi − h. si tomamos xi = 1 y h = 1. Para discretizar las derivadas de una función F (x. demostrar que la fórmula f (xi ) ≈
Problema 52 (2 puntos) Considerar en el problema anterior que xl = xi − h. 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 . y) = (Fx (x. lj). y). y) = Fxx (x.j ∼ = F (hi. 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. utilizando la expresión anterior para aproximar la derivada de f (x) = x3 en x = 1.2 que está más próximo al valor real que utilizando la primera fórmula. y)).
. y el Laplaciano ∆F (x. demostrar que la fórmula f (xi ) ≈ 2
∂ 2 F (x. en este caso. en este apartado. más preciso es el valor de la derivada. si xr = xi + h. y). Problema 49 (3 puntos) Dados 3 puntos distintos xl . F (x − h. supondremos que la dimensión es 2.
5.13 − 0. xi . y comprobar que da la misma fórmula que utilizando los desarrollos de Taylor. Por ejemplo. F (x + h. y) + Fyy (x. calcular la derivada de ese polinomio en xi .
j+1 + Fi−1.j−1 − 4Fi. 0) conociendo los siguientes valores: F (0.j+1 − Fi. consideremos una función tal que en un entorno de un punto (hi0 . F (0. 0) = 0.j−1 + Fi−1. 0) = 0.j−1 ) +γ 4h (Fi. utilizando el desarrollo de Taylor. 1 ) = 1 .j h2
Problema 56 (2 puntos) Calcular una aproximación del laplaciano de una función F (x. F ( 2 . 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.j 2h2 Fi+1. − 1 ) = 1 . 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. F (0. F ( 1 .j+1 + Fi−1.j−1 − 4Fi. hj0 ).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. supondremos que l = h. pueden utilizarse diferentes esquemas.j−1 ) + 2h (Fi+1. de donde deducimos. y). si queremos que ambos valores de ∆F coincidan.j − Fi−1.j + Fi−1. F (− 1 .j+1 + Fi. y) = (0. F (− 1 . 0) = 1 .j ) + 2h (Fi+1. 2 2
.j + Fi−1. lo será en particular para rotaciones de 45 grados. 2 4 2 4 2 4 2 4 F ( 1 . por tanto.j−1 − Fi−1. Para ello. Teniendo en cuenta que la norma euclídea del gradiente es invariante por rotaciones. F ( 1 . 2 ) = −1. un parámetro a elegir.j+1 + Fi−1.j+1 + Fi+1. Hablando en términos 3 de teoría de la señal. Problema 55 (3 puntos) Demostrar. Para simpliﬁcar. y) en el punto (x.j =γ + 2h2 Fi+1. 2 ) = 0. 0) = − 1 . − 1 ) = 1. Por lo tanto.j+1 − Fi−1. hj0 ) = γ 2 + (1 − γ) 2 2h h Ahora bien. es decir: ∆F = Fi+1. 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 . F ( 1 . − 1 ) = 2 2 2 2 2 2 2 1 1 1 1 1 1 1 2 . utilizando el mismo √ argumento que para el ∆F . que γ = 2 − 2. 0) = 1 .j+1 − Fi−1. F (− 2 . que las siguientes expresiones son discretizaciones del laplaciano: ∆F ∆F = = Fi+1. y) = (0.j + Fi. F (− 2 .j−1 − 4Fi.j−1 + Fi+1.Discretización del Laplaciano Para discretizar el operador ∆F en un entorno de 3 × 3 puntos.j )x = (1 − γ) (Fi. F (0. cualquier promediado de las dos expresiones también es una discretización del laplaciano. hj0 ) = γ 2 + (1 − γ) 2 2h h Por lo tanto. − 2 ) = 0. obtenemos la siguiente expresión para el gradiente: (Fi+1.j +(1 − γ) + h2 +O(h) donde γ es un parámetro libre a elegir. 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+1 + Fi+1. 0) conociendo los siguientes valores: F (0. si rotamos 45 grados. 1 ) = 1 . de nuevo.j−1 ) +γ 4h donde γ es. − 1 ) = 1 . F (− 1 . F (− 1 . − 1 ) = 2 2 2 2 2 2 2 2 2 2 2 1 2. 0) = 1 . 1 ) = 1 . 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 . la función inicial en torno al punto (hi0 .j + Fi.j+1 + Fi. 1 ) = − 1 . F ( 1 . Discretización del gradiente Siguiendo el desarrollo de Taylor mostrado anteriormente. Problema 58 (2 puntos) Calcular una aproximación del gradiente de una función F (x. y) en el punto (x.j−1 − 4Fi.j+1 + Fi−1.j+1 − Fi+1. debemos elegir γ = 2 . F (0.
b] utilizando la evaluación de f (x) en ciertos puntos de [a. y para n = 2. Es decir Z b N X P (x)dx = wk P (xk )
Problema 60 (2 puntos) Se consideran.. b] = (−∞. x1 = 0 y x2 = 0. . y dado un intervalo [a. se utilizan los ceros de los denominados polinomios de Hermite..8888888889 − 0. L1 (x) = x.. deﬁnidos como H0 (x) = 1. 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.3478548451 0. 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.6251451549 − 0. Deﬁnición 2 Una fórmula de integración numérica se denomina exacta de orden M si. ∞). hay que emplear otros métodos para aproximar las integrales.6251451549 −0. En este caso. una fórmula de integración numérica se puede escribir como Z
Problema 59 (2 puntos) Aproximar el valor de la siguiente integral.7745966692 0.7745966692 0. Es decir. y Hn (x) = 2xHn−1 (x) − 2(n − 1)Hn−2 (x) para n ≥ 2. para cualquier polinomio P (x) de grado menor o igual que M. a = −∞ o b = ∞. es decir. b] cualquiera.5. los puntos x0 = −0.5 y los pesos w0 = w1 = w2 = 2/3. 1]. 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. −0. 1]. Usar esta fórmula de integración para calcular númericamente la siguiente integral y compararla con el resultado análitico (exacto). la fórmula de integración numérica aproxima la integral de la siguiente forma: Z
. 1]. ¿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. Demostración [Hu] Pg. b]. b].3399810436 0. b].5773502692 1 3 0. la fórmula es exacta.8611363116 0. Z
Deﬁnición 3 Se denominan polinomios de Legendre Ln (x) a la familia de polinomios dada por L0 (x) = 1.Integración Numérica Métodos de Cuadratura de Gauss Sea f (x) una función deﬁnida en un intervalo [a.8611363116 0. cuál sería la fórmula de integración numérica de Legendre utilizando un sólo punto de interpolación.N los ceros del polinomio de x Legendre LN (x). Estos puntos y estos pesos se utilizan para aproximar la integral de una función en [−1. 3..5555555556 4 0. b] es inﬁnito.5773502692 1. 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. En el caso [a.3478548451
Cuando el intervalo [a.3399810436 0. para el intervalo [−1. vamos a aproximar el valor de la integral de f (x) en [a. 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. 0. nLn (x) = (2n − 1)xLn−1 (x) − (n − 1)Ln−2 (x)
Problema 61 (2 puntos) Encontrar. H1 (x) = 2x. utilizando los ceros y pesos asociados a los polinomios de Legendre.5555555556 0...
y Ln (x) = (2n − 1 − x)Ln−1 (x) − (n − 1)2 Ln−2 (x). Una alternativa consiste en dividir previamente la integral en subintegrales de la manera siguiente: Z
Problema 65 (2 puntos) Aproximar. 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. 585 786 438 0. la fórmula de integración numérica aproxima: Z
donde a = x0 < x1 < . para aumentar la precisión es necesario aumentar el grado de los polinomios. se pueden utilizar los desarrollos a partir de los polinomios de Legendre. lo cual resulta complejo para valores grandes de N. 772 453 851 2 −0. para n ≥ 2. A continuación se aproxima numéricamente cada una de las integrales Z xk+1 f (x)dx
Para ello. 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 .
donde a es un número real cualquiera. 1. utilizando dos puntos de aproximación. 853 553 390 3 3. 1.
Fórmulas de Integración Numérica Compuestas Con las fórmulas que hemos visto hasta ahora. 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. 707 106 781 0. 707 106 781 0. ∞). < xM+1 = b. En este caso. el valor de la integral: Z ∞ 1 dx 1 + x2 −∞
Para el intervalo (0. Demostración [Hu] Pg....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 (−∞. 886 226 925 5 0. 2 0. ∞). L1 (x) = 1 − x. se utilizan los polinomios de Laguerre Ln (x). Demostración [Hu] Pg. 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. deﬁnidos por L0 (x) = 1.
utilizando los polinomios de Hermite. ∞). 211-213
. 146 446 609 3
utilizando los polinomios de Laguerre. 414 213 562 0.
la integración a partir de los ceros de los polinomios de Legendre no puede utilizarse.
donde debemos elegir los puntos (xi . 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. Para realizar esta elección se utilizan técnicas de cuadratura. la integral Z 1 ¡ 3 ¢ x − x4 dx
. yj ) y los pesos wij . Para simpliﬁcar la exposición. 2 horas)
Ahora bien. Integración numérica en dimensiones superiores En esta sección. y)dxdy
Por tanto. 2 2
Práctica 5 (Implementación Método de Integración de Simpson. sarrollo en serie de Taylor centrado en el punto xm = xk +xk+1 .j
Aunque estas fórmulas sean menos precisas que las deducidas a partir de los ceros de los polinomios de Legendre. 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. yj )
Ω i. por el método de Simpson. −∞ e−x dx = π = 1. cuando sólo conocemos f (x) en un conjunto de la forma xk = x0 + hk. supondremos que la dimensión es 2. pretendemos aproximar Z F (x. 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 )
Ω i. 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. y)dxdy ≈ wij F (xi . 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. 0 sin(x)dx = 2 Z xk+1 f (x)dx ≈ xk R1 x ¶ Z xk+1 µ 2. − 1 . en este caso. es decir. Es decir. 1 y 1. 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. 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. 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.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 . Nótese que. Es decir: 2 Rπ 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. estudiaremos las técnicas de integración numérica sobre dominios Ω de dimensión superior a 1. 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. 0. La función a integrar se deﬁnirá aparte (como ³ ´ Z xk+1 en el caso del mérodo Müller).
En función de los vértices. 1). y2 ). (1. 0).
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. en este caso. 0) y (0.donde m y n determinan el grado de los polinomios. yk ) x e
F (x. al igual que en dimensión 1. 1]. b]. viene dada. por tanto. el cálculo es un poco más complejo. en general. la exactitud en dimensión 2 la podemos deducir a partir de la exactitud en dimensión 1. Problema 69 (3 puntos) Deducir la fórmula de integración numérica sobre el rectángulo [−1. podemos escribir: Z xm y n dxdy =
AREA(T ) el área del triángulo T. Problema 70 (2 puntos) Deducir la fórmula de integración numérica sobre un rectángulo [a. y [c. d]. En este caso. y1 ). En el caso de que Ω sea un triángulo. el área viene determinada por ¯⎞ ⎛¯ ¯ 1 1 1 ¯ ¯ ¯ 1 AREA(T ) = ABS ⎝¯ x0 x1 x2 ¯⎠ ¯ ¯ 2 ¯ y0 y1 y2 ¯
y. también podemos extender los resultados al caso en que los intervalos sean inﬁnitos. como hemos visto anteriormente. µ ¶ Z ∼ F x0 + x1 + x2 .
. 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. Denotaremos por
Problema 73 (2 puntos) Se considera el triángulo T de vértices (0. b]x[c. 1] resultante de aplicar la integración numérica en una variable en los intervalos [−1. Nótese que. y) = ax + by + c. y) ∼ AREA(T ) =
wk F (ek . y0 )w0
sea exacta para polinomios de grado 1 en x e y. y) ∼ AREA(T ) =
wk F (ek . d]. por los polinomios de Legendre. De estas relaciones se puede deducir. y0 ) y el peso w0 para que la fórmula de integración numérica: Z F (x. d] resultante de aplicar la integración numérica en una variable en los intervalos [a. y) en 4 puntos. y) = 3 3 T Integración sobre triángulos utilizando 3 puntos. y0 + y1 + y2 AREA(T ) F (x. Consideremos un triángulo T de vértices (x0 . y)e−x
F (x. y [−1. y)dxdy ≈ F (x0 . (x1 . 1]. Deducir cual debe ser el punto (x0 . los valores de los puntos y los pesos. x3 e
utilizando integración numérica. b]x[c. yk ) x e
x2 e tos. Z F (x. Es decir P (x. (x2 . 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. que. 1]x[−1. y0 ). Un caso particularmente sencillo es cuando Ω es un rectángulo [a.
3 puntos y 4 puntos. x = (x1 . p = 2. k una norma vectorial. La deﬁnición más utilizada es la denominada norma p.. una norma mide la magnitud o tamaño de un vector x. que corresponde a la norma euclídea. incluyendo técnicas iterativas de resolución de sistemas de ecuaciones y cálculo de autovalores. k x k2 < 1 3. cuando trabajamos en varias dimensiones. esta desigualdad es cierta por la propia deﬁnición de k A k.. ANÁLISIS NUMÉRICO MATRICIAL II En esta sección veremos algunos aspectos más avanzados del análisis matricial. en el espacio vectorial de los números reales. la desigualdad anterior es equivalente a k Ax k ≤k A k kxk Ahora bien. puesto que k x k> 0. utilizando 1 punto. esto es.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 . Problema 79 (2 puntos) Demostrar que si A y B son dos matrices de dimensión N xN. 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. y demostrar que la norma k x kp veriﬁca las propiedades de la deﬁnición de norma. donde p es un número real positivo.
. . 2). entonces. Sin embargo. Si x 6= 0. la desigualdad es trivial. la norma ”natural” es el valor absoluto. deﬁnida por k x k∞ = max | xi |
Problema 75 (4 puntos) Tomar N = 2 . 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. lo que da lugar a la denominada norma inﬁnito.
Normas de vectores y matrices Deﬁnición 4 Una norma k . 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 . se podría deﬁnir su norma considerando la matriz como un vector de dimensión N xN . 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 . y ∈ E.. (2. Básicamente. subordinada a la norma vectorial k . Otro caso interesante es aquél que se produce cuando hacemos tender p hacia inﬁnito. x2 . k una norma vectorial. k x k1 < 1 2. 0).. para cualquier vector x se veriﬁca que k Ax k≤k A k · k x k Demostración: Si x = 0. Empezaremos recordando algunos conceptos relacionados con los autovalores. Entonces. 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 . Por ejemplo. x2 ) que veriﬁcan que 1. Se deﬁne la norma de A. k como k A k= sup
donde Ω es el triángulo de vértices (0. xN ). 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. existen múltiples formas de deﬁnir una norma. Sin embargo. 0) y (0. entonces.
y por tanto el resultado del Teorema. por tanto k Ax k k λx k = = |λ| ≤k A k kxk kxk Lo que demuestra el teorema. 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 . 53. 73. se cumple q k x k2 = (η 1 )2 + .. de la forma: x = η 1 x1 + η2 x2 + . + η N xN Al hacer Ax. además. + (ηN )2 q 2 2 k Ax k2 = (η1 λ1 ) + . 1 e inﬁnito de la matriz µ ¶ 1 0 A= 1 1
. En consecuencia. Sea x un vector cualquiera. Puesto que A posee una base ortonormal de autovectores
Demostración: [La-Th] Pg..Deﬁnición 6 Un autovalor de A es un número λ real o complejo tal que existe un vector x. tendríamos la igualdad. k una norma vectorial.75. Teorema 19 Si los autovectores de una matriz A de dimensión N xN forman una base ortonormal de RN . Teorema 21 Sea A una matriz cualquiera. entonces existe un autovector x tal que Ax = λx. por tanto. Demostración: [La-Th] Pg. + (ηN λN ) k Ax k2 ≤ ρ(A) k x k2 para cualquier vector x. Deﬁnición 8 Se deﬁne el radio espectral de una matriz A como ρ(A) = max{| λi | : λi autovalor de A}
Y. entonces todos sus autovalores son reales y. denominado autovector. 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 (λ). Entonces k A k≥ ρ(A) Demostración: Si λ es un autovalor de A. obtenemos que
Deﬁnición 7 Se denomina polinomio característico P (λ) de la matriz A..
Problema 82 (2 puntos) Calcular las normas 2. sus autovectores forman una base ortonormal de RN . al tomar el supremo en x. + η N λN xN Como los autovectores son ortonormales. tal que Ax = λx
xi . además. el vector x se podrá expresar como una combinación lineal de autovectores. y puesto que los xi son autovectores. la desigualdad se mantiene. en primer lugar. Vamos a demostrar la desigualdad k A k2 ≤ ρ(A) Dado que el teorema anterior determina la desigualdad en el otro sentido.. su producto escalar veriﬁca que (xi . entonces k A k2 = ρ(A) Demostración: Recordamos. al polinomio dado por el determinante P (λ) =| A − λI | que
Ax = η1 λ1 x1 + η2 λ2 x2 + . 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.
Teorema 18 Sea A una matriz y k . 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. xj ) =
Teorema 20 Si una matriz A de dimensión N xN es simétrica.
−1. 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. 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. 1). 0. Problema 84 (2 puntos) Demostrar que. 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. 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. Vamos a considerar ahora el mismo sistema. Problema 85 (2 puntos) Calcular el condicionamiento para la norma 2.84. Es decir. ° ° kδuk ≤ °A−1 ° kδbk Por otro lado. Como podemos observar. 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. de donde δu = A−1 δb y. al mismo tiempo.29. Para ilustrar de qué estamos hablando.2.1). mejor comportamiento numérico tendrá la matriz A.29 = 3029 0.
. −12. 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. multiplicando esta desigualdad con la anteriormente obtenida para kδuk . 3. concluimos la demostración del teorema. se obtiene que Aδu = δb.
lo cual indica un condicionamiento bastante malo.Problema 83 (2 puntos) Demostrar la siguiente igualdad: ρ(t AA) = ρ(A ·t A) Teorema 22 Sea A una matriz cualquiera. para la norma 2.6.9 ⎟ ⎟=⎜ ⎟ ⎠ ⎝ 33.9 ⎛
La solución de este sistema es (9. si los autovectores de una matriz A de dimensión N xN forman una base ortonormal de RN .6. 80. entonces. entonces Limn→∞ k An k= 0 ⇐⇒ ρ(A) < 1 Demostración:[La-Th] Pg. y 30.1 ⎠ 30. también se cumple que kbk = kAuk ≤ kAk kuk de donde obtenemos que 1 kAk ≤ kuk kbk Así. 4.01
cuya solución es (1. 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. 1. 1. cuanto más pequeño sea χ(A). por tanto el condicionamiento sería χ(A) = 30. Obviamente. 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. Consideremos de forma genérica un sistema de ecuaciones de la forma Au = b y.01.5. por tanto.86.1 ⎟ ⎜ 22.
0 −. 707 107 2. 0 ⎟ Rpq (α) = ⎜ 0 . se veriﬁca −1 que (Rpq (α)) =t Rpq (α). Para anular un valor a0 . 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 α ∈ − π .0
De donde α = − π . . si A es una matriz simétrica. dadas dos matrices A y R. cos α .Método de Jacobi Este método se aplica a matrices reales y simétricas. Se basa en el hecho de que. Este método intenta diagonalizar A realizando transformaciones del tipo R−1 AR. 707 107 ⎠ −. En el método de Jacobi se utilizan las denominadas matrices de rotación. . sólo se ven afectadas las ﬁlas y columnas de índices p y q. . Problema 86 (2 puntos) Sean las matrices A y R. y simpliﬁcar el algoritmo. . cot(2α) = cos 2α (aqq − app ) = sin 2α 2apq
Para convertir en 0 el elemento a12 = −1 de la matriz. haciendo 0 los elementos no diagonales mayores en módulo. 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. debemos elegir α tal que cot(2α) = (a22 − a11 ) =0 2a12
t y al hacer la operación R12 AR12 obtenemos ⎛ ⎞ 1. 707 107 ⎝ 0 3. 1 0 ⎠ 0 0 0 0 0 0 1
Es decir. sin α . 0 ⎟ ⎟ ⎜ ⎜ 0 . Concretamente. que son costosas computacionalmente. 1 . 0 ⎟ ⎟ ⎜ ⎝ 0 . π . Al ser una matriz de rotación. la matriz R12 es 4 ⎞ ⎛ √ 1 1 − √2 0 2 1 1 √ 0 ⎠ R12 = ⎝ √
Para evitar tener que evaluar funciones trigonométricas. la matriz A0 también es simétrica. cos α = 1 p 1 + tan2 (α) sin α = tan(α) cos α
donde los cosenos y senos están situados en las columnas y ﬁlas p y q. sign(x) = 1 si x ≥ 0 y sign(x) = −1 4 4 si x < 0. 707 107 −. . Además. cos α . ⎟ ⎜ ⎜ 0 .0 0 −. Al realizar la operación A0 =t Rpq (α)ARpq (α). basta con elegir α tal que pq (app − aqq ) sin 2α + apq cos 2α = 0 2
. Por tanto. q
El método de Jacobi se basa en ir modiﬁcando la matriz A mediante el procedimiento anterior. 0 ⎟ ⎟ ⎜ . se veriﬁca que los autovalores de A son los mismos que los autovalores de R−1 AR. − sin α . 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. . q apj sin α + aqj cos α j 6= p. que tienen la forma siguiente: ⎞ ⎛ 1 0 0 0 0 0 0 ⎜ 0 1 . . 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.
+ T ∗ T ) SI = CO ∗ T PARA j = 1.. N max HACER p=2 q=1 R = ABS(A(p... 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) A(p. p) = A(p. como Ri es una matriz de rotación del tipo Rpq (α)./SQRT (1. Por tanto. .. p) = 0 FIN PARA n PROCEDIMIENTO TERMINADO INCORRECTAMENTE NÚMERO DE ITERACIONES MÁXIMO EXCEDIDO Veamos ahora cómo podemos calcular los autovectores. RM .. j)) p=j q=i FIN IF FIN PARA j FIN PARA i IF R < T OL HACER PROCEDIMIENTO TERMINADO CORRECTAMENTE.. q
Veamos ahora la convergencia del método de Jacobi para el cálculo de autovalores. en cada iteración haremos B = B · Ri . 576-577. Sea A1 = A. · R1 AR1 · .. j) = SI ∗ D + CO ∗ A(q. p))/(2 ∗ A(p. ... Teorema 24 Sea una matriz A simétrica. + C ∗ C)
FIN IF CO = 1. q) − A(p. Los parámetros de entrada son la matriz simétrica A. j) FIN IF FIN PARA j A(p. N
. .. N i = 1.. q) A(q. inicializamos B a la identidad antes de entrar en el bucle. Al utilizar el método de Jacobi.. q)) IF C < 0 HACER T = −C − SQRT (1.. q j 6= p. j) A(j. Además los elementos no diagonales de A convergen hacia 0. p) − T ∗ A(p. de la expresión anterior se obtiene que Abi = dii bi Es decir. entonces los elementos diagonales de la matriz Ak convergen (k → ∞) hacia los autovalores de la matriz A. A continuación. Denotemos por B la matriz B = R1 · . LOS AUTOVALORES SE ENCUENTRAN EN LA DIAGONAL DE A. i − 1 HACER IF ABS(A(i. el número máximo de iteraciones N max y la tolerancia T OL para decidir cuándo son ceros los elementos no diagonales... . DIM HACER IF ( j 6= p AND j 6= q) HACER D = A(p. . q)) PARA i = 3.. j) A(j. para calcular la matriz B en el algoritmo anterior que calcula los autovalores. q) = A(q. · RM = D
donde D es una matriz diagonal que contiene los autovalores de A en la diagonal. j)) > R HACER R = ABS(A(i. q) = A(q.Ahora bien.. añadiremos en cada iteración las operaciones necesarios para ir obteniendo B. q) = A(q..Problema 89 (3 puntos) Dentro del método de Jacobi para el cálculo de autovalores. Algoritmo del Método de Jacobi para el Cálculo de autovalores. de tal forma que
−1 −1 RM · . + C ∗ C) ELSE T = −C + SQRT (1. j) = CO ∗ D − SI ∗ A(q.. q) + T ∗ A(p. PARA n = 1. Numéricamente. y sea Ak la matriz transformada de Ak−1 . su dimensión DIM. Despejando de la anterior igualdad obtenemos que AB = BD Si denotamos por bi el vector columna i de la matriz B...... vamos transformando la matriz A multiplicándola por una secuencia de matrices de rotación R1 . SALIR FIN IF C = (A(q. la matriz B determina los autovectores. Demostración: [La-The] Pg. que se transforman de la siguiente manera: b0 ip b0 iq = cos(α)bip − sin(α)biq = sin(α)bip + cos(α)biq i = 1. haciendo cero el elemento no diagonal mayor en módulo de la matriz Ak−1 . bi es el autovector de A asociado al autovalor dii .. · RM . . 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. En primer lugar.. DIM HACER PARA j = 1.. p) = A(p. ..
0 . utilizando la fórmula :
N 1 X ABS(U (i) − V (i)) ERROR_V ECT ORES = N i=1 ABS(U (i)) + 1.0 −1. ⎝ 1 ⎠ ↔ 2 ⎩ ⎭ 1 ⎛ ⎞ 2 −1 0 2.B. 465 ↔ ⎜ ⎟ ⎜ −. 241 2 ⎞ ⎞ ⎛ . 853 75 ⎟ ⎟ ⎜ ⎜ . Las matrices de dimensión la asignatura.N) : Devuelve la diferencia entre los vectores U y V.0 −2. tomando T OL = 0. 575 36 ⎟ ⎟ ⎜ − 2. 324 6 ⎟ ⎟ ⎜ 1. 06 ↔ ⎜ ⎜ −. 562 8 ⎟ ⎟ ⎜ ⎜ 0.0 −1.TOL. 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. 626 28 .
3. 268 11
⎜ ⎜ ⎜ 16. si deﬁnimos la secuencia un = A se veriﬁca que Limn→∞ sign
¡¡ n n−1 ¢¢ u . 855 2 ⎜ ⎝ 1. 363 27 . el algoritmo de Jacobi se puede modiﬁcar haciendo cero el primer elemento apq que se encuentre que veriﬁque |apq | ≥ T OL. utilizando la función ERROR_VECTORES(). 176 3 1. 165 3 ⎞ ⎛ −1. los vectores A¯i y λi xi .Nota. 826 8 ⎟ ⎟ ⎟ ⎜ ⎠ ⎠ ⎝ 1.AUTOVALORES. 12 ↔ ⎜ −1. B es una matriz donde se guardan los autovectores por columnas. de dimensión N.Niter): Realiza el cálculo de los autovalores y autovectores de una matriz simétrica A por el método de Jacobi.0001 y N iter = 1000: ⎛ ⎞ 2 2 −2 1. 11 ↔ ⎜ ⎟ ⎜ −. 088 8 ⎟ ⎟ 5. ⎝ −1 ⎠ ↔ ⎩ ⎭ ⎩ ⎭ 1 1 ⎞⎫ ⎧⎛ ⎨ 0 ⎬ 4. 796 84 ⎟ ⎟ ⎜ ⎜ 3.N
max ERROR_V ECT ORES(A¯i . A = ⎝ −1 2 −1 ⎠ 0 −1 2 Resultado:
Nota: Obsérvese.Nmax. 942 ↔ ⎜ ⎜ −1. 659 86 ⎛ . 938 54
Comprobar los resultados obtenidos con los siguientes ejemplos. La función devuelve 0 si termina bien y −1 en caso contrario. Para no tener que buscar en cada paso el máximo de los elementos no-diagonales de Ak .Nmax): Para comprobar que los autovalores λi y su autovectores xi están bien estimados. A = ⎝ 2 1 1 ⎠ −2 1 1 Resultado: ⎫ ⎧⎛ ⎧⎛ ⎞ ⎞⎫ 1 ⎬ ⎨ −2 ⎬ ⎨ ⎝ −1 ⎠ ↔ −2. 984 4 ⎜ −. 901 1. 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. 985 4 ⎟ ⎟ ⎜ ⎠ ⎝ 1. λ¯i ) x x
⎞ −.05 170 ⎟ ⎟ ⎜ 2. ⎝ 2 ⎠ ↔ 2 − 2 ⎭ ⎩ 1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎛
• FUNCTION JACOBI(A.0 . 741 24 ⎜ ⎜ −. comparar ¯ para cada autovalor λi .N.
⎞ 2. entonces.0 . ⎛ 0 1 6 0 0 0 ⎜ 1 0 2 7 0 0 ⎜ ⎜ 6 2 0 3 8 0 4.u k un k= λmax
.N. ⎝ − 2 ⎠ ⎭ ⎩ ⎭ ⎩ 1 1 ⎧⎛ ⎞⎫ 1 ⎬ √ √ ⎨ √ 2. 746 5 ⎜ .AUTOVECTORES. • FUNCTION ERROR_AUTOVECTORES (A. A = ⎜ ⎜ 0 7 3 0 4 9 ⎜ ⎝ 0 0 8 4 0 5 0 0 0 9 5 0 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠
⎧⎛ ⎧⎛ ⎞⎫ ⎞⎫ 1 ⎬ ⎨ −1 ⎬ ⎨ √ ⎝ 0 ⎠ ↔ 2. 896 18 ⎟ ⎟ ⎟ ⎜ ⎟ ⎟ ⎜ ⎟ − 12. al comparar los resultados. Sea un vector u1 no ortogonal al subespacio engendrado por los autovectores del autovalor λmax . 215 57 ⎟ ⎠ ⎝ 1.V. que los autovectores están deﬁnidos módulo la multiplicación por una constante. 729 18 ⎟ ⎜ −. 6 ↔ ⎜ ⎜ ⎜ ⎝ ⎛
. 593 8 ⎜ −1. 521 32 ⎜ . 400 89 ⎟ ⎠ ⎝ 1. 998 09 −10. x ¯ Devolver la expresión ERROR_AU T OV ECT ORES =
Sean x1 . si λmax es negun = A−1 k un−1 k ativo. veriﬁca que Además ¡¡ ¢¢ 1 Limn→∞ sign un .. dicho autovector tiene norma 1. el término sign un .. En este caso.u k un k= λmax ¡ ¡¡ ¢¢¢n Por otro lado.. µ ¶ 2 1 A = 0 1 µ ¶ −3 0 A = 1 1
An−1 u1 µ λn−1 x1 + . entonces u1 se puede escribir como u1 = µ1 x1 + . Por tanto. En primer lugar.. un−1
Además. si λmax es positivo.u ≥ 0 y sign un .¡¡ decir sign un .¡ ¡¡ ¢¢¢n Limn→∞ sign un . todos los co1 cientes de la forma λmin = 0 max{λi autovalores de A−1 } µ ¶n λi Por tanto. un−1 = un = k un k
Para n = 1 la igualdad se cumple por la deﬁnición de u2 . |λmax | obtenemos que la secuencia tienden hacia 0.. y demostrémoslo para n: un =A k un k
µ1 x1 + . salvo si λi = λmax . un−1 = −1 si un . un−1 = 1 si es ¡ n n−1 ¢ ¢¢ ¡ ¢ u ... dicho un−1 cociente es 1n . + µN |λλN | xN ° |λmax | ° ° max = |λmax | y por tanto ¡¡ n n−1 ¢¢ u . + µN xN donde supondremos que x1 es un autovector asociado a λmax y que µ1 6= 0.. + µM xM kµ1 x1 + . un−1 para n suﬁcientemente grande es 1n si λmax es positivo o (−1)n si λmax es negativo.. + µN |λmax | xN ° |λmax | ° ° °= = |λmax | ° ³ ´n−2 ³ ´n−2 ° ° °µ1 λmax x1 + . un−1 < 0. xM los autovectores asociados a λmax . + µN |λλN | xN ° culo del autovalor de módulo menor λ . u )... un−1 es el signo del¢¢ producto ¡¡ escalar de un y un−1 . obtenemos que Limn→∞ sign ¡ ¡¡ ¢¢¢n Limn→∞ sign un ... Supongamos que se cumple para n−1... + µ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 + . Por la igualdad anteriormente demostrada obtenemos que
que es un autovector de λmax de norma 1. 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 + .. o (−1)n . Demostración. 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.. ¡¡ ¢¢ Teorema 26 sign un .. ... + µM xM k
Con lo que queda demostrado este primer resultado.. como A posee una base de autovectores. si aplicamos el método anterior a A−1 . + µN λn−1 xN max N °= = ° 1 n−2 °µ1 λmax x1 + .. realizando 3 iteraciones en el método. para n suﬁcientemente grande el signo de n n−1 λmax coincide con el signo del producto escalar (u . Por otro lado. que denotaremos por xi . teniendo en |λmax | ° ° max min cuenta que Cuando hacemos tender n hacia inﬁnito. 1). y u1 no es ortogonal al espacio generado por los autovectores asociados a λmax . hasta calcular u4 y partiendo de u1 = (1. un−1 k un k= λmin un =
podemos perturbar ligeramente el valor de µ para que IGAUSS() no dé problemas. Si hacemos iteraciones del esquema anterior a partir de la aproximación inicial u1 = (0. 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. Por ejemplo. entonces se obtiene que el autovalor menor de A0 es justamente λ. y. Para autovalores que se encuentren entre λmin y λmax . haremos 1 = 10−11 A = A − µId J = IGAU SS(A0 .. .. quedaría como sigue: Si µ es el autovalor que estamos tratando. donde µ es uno de los elementos diagonales de la matriz que resulta de aplicar el método de Jacobi. por tanto. entonces el autovalor más pequeño de A0 está muy próximo a 0. 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. obtenemos que x2 y2 z2 = = = 1+0 1 = 2 2 0+0 =0 2 1+0 1 = 2 2
. Nótese que si el autovalor µ está calculado con mucha precisión. Algorítmicamente. 0. 1).. se evita calcular directamente A−1 . 1).Limn→∞ sign
un es un autovector de λmin k un k En los casos prácticos. utilizando el método de Jacobi.. y1 .h IGAUSS(). al hacer iteraciones de la forma un = M un−1 + c se obtenga que un converge hacia u.. y como el determinante de una matriz es el producto de sus autovalores.0) T HEN = ∗ 10. 1. 1. µ = µ(1 + ) GOT O 1 EN DIF
Buscar la solución de este sistema es equivalente a buscar un vector u = (x. Para evitar esto. 1). 1) y tomando como norma kuk = maxi |ui |. si consideramos la matriz A0 = A − µI. 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. la solución del sistema original.
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. y. la solución exacta del sistema es u = (1.u
Para ello.) IF (J. y se obtiene un resolviendo el sistema Aun = un−1 k un−1 k
¡¡ n n−1 ¢¢ u . 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 . 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. 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. 0). realizando 2 iteraciones del método a partir de u1 = (1. podemos aplicar el método de la potencia inversa anterior.N E. calcular dos iteraciones del método de la potencia inversa partiendo de u1 = (1.
y vamos actualizando las componentes del vector aproximación según las vamos calculando. Todas se basan en descomponer A de la forma A = L + D + U.**120 si P termina no correctamente..96 ⎠ 0... Nf la dimensión real. 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. 1.. básicamente. − a2N un−1 + b2 1 3 N = a22 .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).Nf )
... − a2N un−1 + b2 1 3 N = a22 ...5 . −aN1 un − aN 2 un .Tolf ) que devuelve el autovalor máximo de una matriz y su autovector por el método de la potencia. En este caso. −aN1 un−1 − aN 2 un−1 . − a1N un−1 + b1 2 N a11 −a21 un − a23 un−1 . y U es la matriz triangular superior que corresponde a la parte de A situada por encima de la diagonal. Esta función devuelve el valor 2. Este método consiste en tomar MGS cGS = (D + L) = (D + L)
A efectos prácticos. 1). la aplicación de este método no requiere el cálculo directo de la matriz inversa (D + L)−1 . 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. obtenemos el método de Gauss-Seidel. u8 = ⎝ .uf. y la función AUTOVALOR_MAXIMO(Af. obtenemos que ⎛ ⎞ ⎛ ⎞ ⎛ ⎞ 0. puesto que el paso de una iteración a otra puede hacerse de la siguiente forma: un 1 −a12 un−1 − . después de haber calculado todas las componentes por separado.73 ⎠ . Ejemplo 18 Vamos a aplicar el método de Gauss-Seidel al sistema del ejemplo anterior.84 .Como puede observarse. Por tanto. Tomar como norma kuk = i ABS(ui ) (28 líneas de código como máximo). − 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... las sucesivas iteraciones se van aproximando a la solución u = (1..25 ⎠ . y Tolf la tolerancia. 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. 0. es decir 2x − y −x + 2y − z −y + 2z
Problema 95 (3 puntos) en Fortran las funciones SIGNO_PRODUCTO_ESCALAR(uf. Los parámetros son la matriz Af.Nfmax.Nﬁter. y en el caso de Jacobi se actualiza sólo al ﬁnal...84 . u17 = ⎝ . Nﬁter número máximo de iteraciones... donde D es la matriz diagonal que corresponde a la parte diagonal de A. entonces u veriﬁca que u = Mu + c
De la misma forma.. 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 − . − aNN −1 un −1 + bN 1 2 N = aN N =
Es el que se ha utilizado en el ejemplo anterior.Nf. el vector candidato inicial uf. Nfmax.vf.98 u3 = ⎝ 0..5 .. L es la matriz triangular inferior que corresponde a la parte de A situada por debajo de la diagonal. − a1N un−1 + b1 2 N a11 −a21 un−1 − a23 un−1 . 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. la dimensión para coger memoria.
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. u8 = ⎝ .. y −1 c = (D + U ) b. 0. el siguiente resultado muestra la forma de calcular el valor óptimo de w. −aN1 un .
1 En este caso. u3 = ⎝ .. Se toman. 686 .. Demostración [La-Th]. 0). Ejemplo 19 Vamos aplicar el método de relajación al sistema del ejemplo anterior. 25 ⎠ .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... 0)
Problema 99 (1 punto) Una variante del método de Gauss-Seidel consiste en tomar M = (D + U )−1 (−L). 0. ρ(MJ ) = √2 y wopt = 1. − a2N uN + b2 = w + (1 − w)un−1 2 a22 . Pg.. 785 . 921 . 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. en este caso. − 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. 999 ⎠ . matrices con todos los elementos nulos salvo la diagonal principal y sus codiagonales.5 . de la forma siguiente: un 1 un 2 un N −a12 un−1 − .. un problema difícil.. 0) y tomando w = wopt = 1. 585 . 17. 0. obtenemos que ⎞ ⎛ ⎞ ⎛ ⎞ . el valor de wopt se encuentra siempre entre 1 y 2.. Mw cw = (D + wL)
Si hacemos iteraciones del esquema anterior a partir de la aproximación inicial u1 = (0. 342 ⎠ . en general. en este caso. en el caso de matrices tridiagonales.. obtenemos que ⎛ ⎞ ⎛ ⎞ . Indicar.17. − a1N un−1 + b1 2 N + (1 − w)un−1 1 a11 n−1 n −a21 u1 . qué diferencias de implementación habría con respecto al caso anterior. Sin embargo.. es decir. u8 = ⎝ . 625 . 802 ⎠ . 999 ⎛
. 976 56 u3 = ⎝ .. 976 56 ⎠ . 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..358-362. 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. 988 28
Como puede observarse de la expresión anterior. 999 u2 = ⎝ ..
| aii |>
. Demostración.
| aij | ∀j. por tanto. si w ∈ / (0. 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.
entonces el método de Jacobi asociado al sistema Au = b converge para cualquier aproximación inicial. Teorema 29 Se considera el esquema iterativo un = M un−1 + c.N − aN. ⎜ aN −1. 2). Demostración: En primer lugar.partiendo de u1 = (0. [La-The] Pg.N−1 −1.N −1 aN. 0). − a2N ⎜ a22 ⎜ . Demostración: En primer lugar. 0 − aNN−1.N ⎝ aN−1. − N.1 − aN−1. entonces |1 − w| ≥ 1 y. esto es ⎛ ⎞ 2 −1 0 ⎝ −1 2 −1 ⎠ 0 −1 2 satisface las hipótesis del Teorema anterior. obtenemos que. Por tanto. Entonces en = M n−1 e1 Demostración: La solución del sistema satisface que u = M u + c. como el determinante de una matriz es el producto de sus autovalores. observamos que la matriz MJ puede expresarse como: ⎛ 0 − a12 − a13 .N . obtenemos que
| aij | ∀j. que kMJ k < 1 para la norma 1 o inﬁnito.N−1 . . Ejemplo 20 La matriz del sistema ejemplo tratado anteriormente. − a1N a11 a11 a11 a21 a23 ⎜ − a22 0 − a22 . 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. entonces los métodos iterativos convergen. por las condiciones del teorema.
con la desigualdad estricta en al menos una ﬁla o columna.346-347. un . en consecuencia. . 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.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. Restando esta igualdad de la igualdad un = M un−1 + c. Mw posee al menos un autovalor de módulo mayor o igual que uno.2 a ⎜ − aN −1. Teorema 33 Si A es una matriz irreducible y se veriﬁca que X | aii |≥ | aij | ∀i. 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.1 aN. Además. su determinante es el producto de los elementos diagonales. / entonces ρ(Mw ) ≥ 1. 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. el teorema se concluye teniendo en cuenta que cualquier norma de una matriz es siempre mayor o igual que su radio espectral. 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.2 a − aN. se tiene.
Por lo tanto. . 0. Calcular previamente el parámetro de relajación óptimo. Teorema 32 Si una matriz A veriﬁca que X | aij | ∀i. .N −1 0 aN. observamos que las matrices D +Lw y (1−w)D −wU son matrices triangulares y. 2).
⎝ −1 2 0 ⎠ ⎝ y ⎠ = ⎝ 3 ⎠ 0 −1 3 z 1 ⎛ ⎛ ⎞⎛ ⎞ ⎛ ⎞
donde u0 es una aproximación de la solución de f (u) = 0.. un sistema no lineal de ecuaciones de dimensión N. la dimensión para coger memoria. que por defecto se tomará 0. Si el método no converge devuelve −1.. si una matriz A veriﬁca que por ﬁlas o columnas su suma es siempre igual a 0. uN ) = 0 .Nﬁter) que devuelve el condicionamiento de una matriz utilizando el método de Jacobi para calcular los autovalores.
4. uN ) = 0 donde f (u) = (f1 (u). En general. 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...... el parámetro de relajación w. ⎝ −1 2 −1 ⎠ ⎝ y ⎠ = ⎝ 0 ⎠ 0 −1 2 z 1 ⎞⎛ ⎞ ⎛ ⎞ 1 3 3 x 7 3. Probar el método para los sistemas 1 −1 0 x −1 1.. ⎝ 3 1 3 ⎠ ⎝ y ⎠ = ⎝ 7 ⎠ 3 3 1 z 7 ⎛
. Se supondrá implementada la función JACOBI(A. Los sistemas ejemplos del directorio de la asignatura.. de un polinomio de grado 2 dado por P2 (z) = az 2 + bz + c.. Por ejemplo.Niter) que devuelve 0 si termina bien y 1 si termina mal.TOLf. 2 horas) Desarrollar una función en Fortran 77 donde se implemente el método de relajación.. si u1 y c son combinaciones lineales de autovectores de M correspondientes a autovalores de módulo menor que 1.. un vector u donde se almacenará la solución. Nﬁter número máximo de iteraciones. .N. reales o complejas.TOL. 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. . Estos ejemplos tienen siempre como solución el vector (1. La función CONDICIONAMIENTO devuelve 2... . el número máximo de iteraciones N max. . Los parámetros son la matriz Af. se escribe como N ecuaciones del tipo f1 (u1 .. Los parámetros de la función serán: la matriz A. . f2 (u)..Nf. es equivalente a resolver el sistema ax2 + bx − ay 2 + c = 0 2ayx + by = 0 que es un sistema no lineal de ecuaciones.*120 si termina mal porque Jacobi da un error o se produce una división por cero.Nfmax. y la tolerancia T OL para evaluar la diferencia entre un y un−1 . uN ).
Práctica 7 (Método de relajación.
Método de Newton-Raphson para sistemas de ecuaciones no lineales En las aplicaciones reales. y u = (u1 . 1.. donde z = x + yi.. Problema 102 (2 puntos) Demostrar que. entonces el método converge. aunque el radio espectral de M sea mayor que 1.Problema 101 (2 puntos) Escribir en Fortran la función siguiente: CONDICIONAMIENTO(Af. y por tanto el sistema asociado a A no tiene solución.. y que inicialmente será el vector aproximación inicial. muchas veces nos encontramos con sistemas no lineales de ecuaciones. 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). el vector b. uN ) = 0 f2 (u1 . y Tolf la tolerancia (21 líneas de instrucciones como máximo).Nmax. fN (u)) es una función de <N → <N .. . entonces el determinante de A es cero. La función devolverá el número de iteraciones necesarias para alcanzar la solución. Nf la dimensión real.. Nfmax. calcular las raíces. 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.. fN (u1 . Comparar la diferencia en la velocidad de convergencia entre el método de Gauss-Seidel y el Método de relajación..
Por motivos de coordinación entre los programas teórico y práctico de la asignatura.. Un ejemplo clásico de ello es el desarrollo de Taylor de una función en un punto a. los polinomios de grado N pueden tener un carácter fuertemente oscilante.N . 1). f 0 (a) f N ) (a) (x − a) + .. Interpolación por splines cúbicos Uno de los problemas básicos del polinomio interpolador de Lagrange. 1). 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. Por tanto. para valores grandes de N. que son polinomios de grado menor o igual que (N + 1)(M + 1) − 1 dados por las siguientes condiciones: ½ ∂ l Hi. x1 = 1. sino también el valor de sus derivadas.. 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. − 4 . Por tanto. resulta de interés interpolar no sólo el valor de la función en ciertos puntos {xi }i=0. como indica el ejemplo siguiente. + (x − a)N 1! N!
que ya es una buena aproximación de la solución exacta dada por el vector (0. siendo ésta la segunda parte.
Interpolación de Hermite En ocasiones.. y. 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. u3 viene dado por 40 µ 1 ¶ µ 13 ¶ µ ¶ 3 − 40 − 40 3 4 u = + = 3 9 39
Tomemos como aproximación inicial u1 = (1.j (x). aproximamos f (x) por un polinomio de grado N .Si partimos de u1 = (1.j (x)
A partir de u1 = (1. para obtener u2 tenemos que resolver µ ¶µ ¶ µ ¶ 2 −2 z1 −1 = 2 2 −2 z2 ¡ 3 1¢ que tiene por solución − 4 .N . u2 viene dado por µ ¶ µ 3 ¶ µ 1 ¶ 1 −4 2 4 u = + = 3 1 −1 4 4 Para calcular u3 . Problema 105 (2 puntos) Plantear el algoritmo necesario para calcular. el polinomio interpolador de Hermite se deﬁne como: P (x) =
(xi )Hi. las raíces complejas o reales de un polinomio de grado 3. 1). es que. En este caso. En el caso general.
donde ξ es un valor intermedio entre x y a. 1). y los resultados obtenidos por la interpolación pueden no ser muy satisfactorios. 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.
Problema 108 (3 puntos) Calcular los polinomios base de Hermite que corresponden a tomar como puntos de interpolación x0 = −1.. 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 . el tema de interpolación de funciones se dividió en dos partes. y) = (1. 1). 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. z) = (1. utilizando el método de NewtonRaphson.
. 1. − 40 . se utilizan los denominados polinomios base de Hermite Hi. 1).. calcular u2 y u3 utilizando el método de Newton-Raphson para aproximar un cero del sistema no lineal.. y el orden de derivación M = 1.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.
hay que añadir una ecuación que involucre a c0 y otra ecuación que involucre a cN . −2. cuando se trabaja con muchos puntos de interpolación. cN ) y N −1 ecuaciones.75
i ∂ 2 P3 ∂x2 (xi+1 ). .. c0 ..5
2ci+1 = 6di hi + 2ci de donde. Nótese que.. . 4. −1. N ci+1 − ci = i = 0... 1. Ejemplo 23 (x2 − 1)(x2 − 4)(x2 − 9)(x2 − 16)(x2 − 25) = −14400
las condiciones anteriores. −4.. 3. 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.... Dicho sistema tiene N +1 incognitas (c0 .... tenemos que buscar 4N valores. se imponen las siguientes condiciones:
i P3 (xi ) = f (xi ) i = 0...
2.. lo que signiﬁca que c0 cN
. para deﬁnir los polinomios.5 x 5
1. . xi+1 ]. entonces ai di bi = f (xi ) i = 0. es decir: a0 .. De la condición
3... N − 1 i P3 (xi+1 ) = f (xi+1 ) i = 0. el número de polinomios es N. El primero consiste simplemente en ﬁjar c0 = cN = 0. . que son polinomios de grado 3. 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 − 1.
bi = 3di−1 h2 + 2ci−1 hi−1 + bi−1 i−1 y.. . Para completar dicho sistema. . . Para añadir estas dos ecuaciones hay dos procedimientos estándares. 2. obtenemos que
0 2. d0 . N − 2 ∂x2
= Vamos a introducir la notación hi xi+1 − xi . b0 .... . . .
i Demostración De la condición P3 (xi ) = f (xi ). se suele interpolar la función utilizando polinomios a trozos. . 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.. −3. como veremos posteriormente... 5]. N − 2 ∂x i+1 ∂ 2 P3 (xi+1 ) 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 3hi ai+1 − ai hi (2ci + ci+1 ) = − hi 3
i = 0. de la condición obtiene que
i−1 ∂P3 ∂x (xi )..5
i De la Condición P3 (xi+1 ) = f (xi+1 ). despejando. .. i = 0. 0.. satisface
Nótese que esta última relación determina un sistema de ecuaciones donde las incógnitas son las variables ci . Para deﬁnir estos polinomios... bN−1 . . 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 .. aN −1 . dN−1. Por razones técnicas. xi+1 ]..25 0 -5 -2.. obtenemos que bi = ai+1 − ai − di h2 − ci hi = i hi ai+1 − ai hi (2ci + ci+1 ) − hi 3
Finalmente.. Si hay N +1 puntos... se obtiene de forma inmediata que ai = f (xi ). N − 1
Despejando. vamos a utilizar también los valores aN y cN . cN−1.
i Teorema 34 Si una familia de polinomios P3 (x) = di (x− 3 2 xi ) + ci (x − xi ) + bi (x − xi ) + ai . N. deﬁniendo un polinomio distinto para cada intervalo [xi .Ejemplo 22 El polinomio base de Lagrange centrado en 0 sobre los puntos xi = −5. despejando todo en función de ci . La técnica más conocida es la interpolación por splines cúbicos..
6 + 0 = 2− = 0. al trazar la curva. y [2.25 2. [1. en primer lugar.133 3 = 1−
Por tanto. no es posible distinguir geométricamente.733x3 + 1.5 0. siguiendo con el resultado del teorema anterior.El segundo procedimiento se utiliza cuando utilizamos los valores de f 0 (a) y f 0 (b).933 3
y 5 4 3 2 1 0 0 -1 x -2 -3 -4 0.
2.8 c2 Los valores bi y di vienen dados por d0 d1 d2 = = = −2. En este caso hi = 1.733 3 2. y P2 (x). 2].5 2.75 0. sabiendo que f (0) = 0. 3] : Como puede observarse.5
0.5 2. Veamos ahora gráﬁcamente el perﬁl de la derivada de los polinomios P0 (x).2 = 1.75 2. por las condiciones sobre las derivadas que hemos impuesto.5 1.2 = 1.25 0. En este caso.133 (x − 2)
Por lo tanto. sobre la gráﬁca de la derivada segunda los puntos de unión se detectan en los lugares donde encontramos un pico. 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.467 (x − 1) + 1
P0 (x) = −0.733x
P2 (x) = −0.933 (x − 2) + 2. Es decir.25 2. Debemos deﬁnir 3 polinomios distintos que corresponden a los intervalos [0.75 -1 -1.25 0 -0. P1 (x).75 1
1. 2. Ejemplo 24 Vamos a calcular los polinomios interpoladores utilizando splines cúbicos al interpolar la función f (x) en los puntos x = 0. tampoco sobre la derivada se aprecian los puntos de unión de los polinomios.467 3 5. y 3. imponemos que
2. f (2) = 0.25 0. 1.8 = −1 − = −0.75 1 1. para calcular los splines cúbicos es necesario.75 3
. 1].2 − 0 = −0.75 2 2. Los bj y dj se calculan directamente a partir de las relaciones mostradas en el teorema anterior.4 + 2.8 = −0.5 1.25 -0. tomar ai = f (xi ). Sin embargo. 2. tal y como se muestra en la gráﬁca de la derivada segunda siguiente:
cuya solución es ¶ µ ¶ µ c1 −2. a simple vista. cuales son los puntos de unión entre los tres polinomios. f (3) = 2.2 (x − 1)2 −0.8 + 2. − 2.25 1.733 3 −4. f (1) = 1.2 = . A continuación.25 2 1.75 3 x
como puede observarse.25
0.25 1 0. se resuelve un sistema de ecuaciones tridiagonal para el cálculo de los ci . 3].8 (x − 2) + 0.5 -0.5
1. tomando c0 = c3 = 0.25 1.667 3 0.667 (x − 1)3 − 2.75 1. parece. el trazado de una única función.5 0. los polinomios son P1 (x) = 1.
2 1 0. f (3) = 2. 1 y 2.. la diferencia entre los diferentes tipos de interpolación también lo hace.4 0.5 0.8 1.5 2 2. f (1) = 1.75 3 x
Ejemplo 26 Consideremos la función f (x). interviene
1. la interpolación por splines cúbicos es la menos oscilante. tal que f (0) = 0. utilizando la función sin c(x). la interpolación por splines cúbicos y la interpolación a través de la función seno cardinal.75 1 1.5
0. 2 2 La interpolación trigonométricos a través de polinomios
. de 0.y
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.8 0.5
0 0 0.2 0 0 0.6
Esta función tiene la propiedad de que en x sin c(0) = 1.
0. Dada una función f (x). f (2) = 0. − π .75
0. la función f (x) = sin(x) en los puntos x = −π. π . Por otro lado. y 3. los polinomios de la interpolación por splines cúbicos (línea sólida). deﬁnida en los puntos x = 0.25 1.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. .. 0. 1.5 1. 0. su función polante en los puntos xi = a · i para i = M. 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.25 0 -50 -25 0 25 x 50
2 1. deﬁnida por sin c(x) = cuya gráﬁca es sin(x) x
1. 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). cuando el número de puntos de interpolación aumenta. 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. π.5 2. y la interpolación utilizando la función sin c(x) (línea a trozos). N dada por la función ¢ ¡ sin(π x − i ) e ¡x a ¢ f (xi ) f (x) = π a −i i=M
= 0.25 2.75 2 2. Problema 110 (2 puntos) Calcular la función que interpola.5 1 1.5
0.25 0..6 0.4 1. y para cualquier entero i distinto sin c(πi) = 0.
. . observamos que... tomando N = 2... cN ) = cl e f (x) − eikx dx = 0 ∂ck −π
Problema 111 (3 puntos) Calcular el polinomio trigonométrico. π]..N . en general complejos.. Dado un conjunto de valores {(xi .. la aproximación mínimo cuadrática lineal consiste en buscar la recta y = ax + b.5 x
Demostración En primer lugar. .25 0 1. Aproximación por mínimos cuadrados La aproximación mínimo cuadrática aproxima. 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 ∈ [− π . y por tanto ! Z π Ã N X ∂E ilx (c−N . dada la forma cuadrática del funcional E(c−N .. . Por otro lado.75
0. b) = sea mínima.Dada una función f (x). π]. tal que la función de error cuadrático E(a. cN ) =
0. π ] / 2 2
Vamos a calcular el polinomio trigonométrico interpolante para N = 3. las derivadas parciales de E(c−N .25 0 -2. deﬁnida en el intervalo [−π.. . yi )}i=1... a través de una función. 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 . un conjunto de valores de forma global.5 -1. 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. el polinomio trigonométrico interpolador es
Demostración En primer lugar. que interpola la función f (x) = |x| en el intervalo [−π. observamos que. en un mínimo. cN ). cN ) con respecto a cualquier ck son cero. debe poseer un
.25 2.5
0... 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.. dada la forma cuadrática que tiene el funcional. éste debe poseer mínimos. π ] 2 2 f (x) = 0 si x ∈ [− π . pretendemos aproximar f (x) como f (x) ≈
donde ck son coeﬁcientes. sin exigir que la función aproximante pase exactamente por ese conjunto de puntos.
F.. da una visión general sobre los últimos avances en Análisis Numérico. destaca por una exposición simple y al mismo tiempo clara. 1988. muy reciente. 1993. [Ci] Ciarlet P. b). [Hi] Higham N. John Wiley and Sons. ”Numerical Methods for Engineers and Computer Scientists”. ”Programación en fortran 77” Anaya. y por tanto X ∂E (axi + b − yi ) xi = 0 (a. sin pretender ser exhaustiva. así como el cálculo de autovalores y vectores propios.W. Vol. y cuya resolución lleva al resultado establecido en el teorema. Trae una buena selección de problemas. 2 Méthodes itératives ". 1996. b) = 2 ∂a i=1
X ∂E (axi + b − yi ) = 0 (a. Esta obra. Keller H. 1 Méthodes directes y Vol. Addison-Wesley Iberoamericana. [St] Stewart G. [Ki-Ch] Kincaid D.E. sin pretender ser tan exhaustiva como otras obras de carácter más general. Grupo Editorial Iberoamérica 1985. 1990. En esta obra se presenta el lenguaje de programación fortran 77 con numerosas aplicaciones al análisis numérico. "Análisis Numérico". Esta obra. 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. Además. la buena presentación de los temas elegidos la hacen de interés. Con un exquisito rigor se abordan los temas básicos del Análisis Numérico Matricial y métodos de optimización. 1989. Contiene todos los tópicos habituales con una descripción muy completa y detallada. [La-Th] Lascaux P. Uno de los libros clásicos más conocidos en Análisis Numérico. ”Introduction à l’analyse numérique matricielle et à l’optimisation”. Masson. The Benjamin/Cummings Publishing Company. SIAM. Excelente libro de base para un curso de Métodos Numéricos. incluyendo la resolución de sistemas a través de métodos directos.E. con múltiples ejemplos y una descripción de los algoritmos bien diseñada. [Hu] Hultquist P. "Accuracy and Stability of Numerical Algorithms".. Masson . 1994.
. b) = 2 ∂b i=1
Esto da lugar a un sistema lineal de ecuaciones cuyas incógnitas son a y b. ”Analysis of Numerical Methods”. Cheney W. ”Afternotes on Numerical Analysis” SIAM.E. trata en profundidad todos los tópicos relacionados con el Análisis Numérico Matricial.. 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. Faires D. ”Análisis Numérico”.. Esta obra. Destaca por el rigor matemático en su exposición. haciendo especial énfasis en la precisión de los algoritmos numéricos y en la propagación de errores. Théodor R. [Bo] Borse G. en un mínimo del funcional E(a. dividida en dos volúmenes. "Analyse numérique matricielle appliquée à l’art de l’ingénieur. Los algoritmos están muy bien descritos a través de un seudocódigo. Inc. Esta obra es un clásico del Cálculo Numérico. las derivadas parciales son cero. Su mayor virtud es el rigor matemático con el que se tratan los temas y una cuidada presentación.
BIBLIOGRAFÍA BÁSICA [Bu-Fa] Burden R. 1996 Esta obra. iterativos y métodos tipo gradiente..mínimo. muestra las últimas tendencias en cuanto a la enseñanza de los conceptos básicos del Análisis Numérico.G. 1966. presenta una cuidada selección de temas básicos en Análisis Numérico.
. escritura y ejecución de un ﬁchero. y el modo edición.f practica1. No poner EN D al ﬁnal del programa principal o de la función. un comentario y un ejemplo. aparece el comando UNIX.name :q! sale del vi sin guardar cambios. Hacer > man chmod para mirar las opciones.name escribe el ﬁchero actual en el ﬁchero f ichero. Finalmente. cd (cd) cambia el directorio activo >cd /users/p701/fortran77 more (type) visualiza el contenido de un ﬁchero >more /users/p701/fortran77/programas/prog1. En primer lugar. chown cambia el propietario de un ﬁchero.APÉNDICE A: Resumen de los comandos de UNIX En este breve resumen seguiremos el siguiente esquema. que es donde se escribe normalmente el texto. A continuación. 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. Hacer > man chown para mirar las opciones. entre paréntesis. >cp /users/p701/fortran77/programas/prog1. 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.f .
APÉNDICE C: Algunos fallos comunes en Fortran 1.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. Puede ejecutarse en dos modos. 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. >mv prog1.f ls (dir) visualiza contenido de un directorio >ls /users/p701/fortran77 cp (copy) copia un ﬁchero en otro.name edita el ﬁchero f ichero. rm (del) borra un ﬁchero >del prog1. :w f ichero. Este comando es de utilidad para salvaguardar la información de directorios y ﬁcheros de miradas ajenas. El modo comando (el que está por defecto al entrar en vi).
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. su equivalente en MS-DOS (si existe). donde se ejecutan comandos.f chmod cambia los permisos de lectura.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.
declarar los tipos de todas las variables que se utilicen al principio del programa o función. Solución: Escribir 1. 10. 6. 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.2. Solución: Poner una declaración de PARAMETER al principio de la función y con ella asignar las memorias de forma estática. Las sentencias GOT O pueden diﬁcultar el seguimiento del ﬂujo del programa y sólo hay que utilizarlas cuando sean indispensables. Por ejemplo. los programas pueden fallar por errores de redondeo en los cálculos.
. Utilizar variables enteras como ﬂotantes o al revés. si hacemos ”f77 -rn prueba. 4.f -o prueba” donde n es 8 ó 16./2. 7. No respetar los tipos en los pasos de parámetros de las funciones. Fortran da la posibilidad de cambiar el número de bits utilizados para almacenar las variables en el momento de la compilacion. No pasar la dimensión de un vector como parámetro de una función. 9. Exceso de sentencias GOT O. No poner ningún comentario en los programas. Utilizar vectores sin declararlos con la sentencia DIMENSION. 8. Siempre hay que buscar que el número de anidamientos sea mínimo. aumentaremos la precisión de la aritmética para las variables reales. Escribir números como 1/2 ó 10 ∗ ∗20 en precisión entera. 11. ∗ ∗20.ó 10. si en lugar de utilizar la directiva -rn utilizamos -dn aumentaremos la precisión de las variables declaradas DOUBLE PRECISION. 5. Análogamente. Anidar excesivamente los programas. A veces. y si utilizamos -in las variables enteras. Sugerencia: Aunque no sea necesario. 3.
Documents similar to Apuntes+de+analisis+numericoMÉTODO DE MULLERSolucion de Sistemas de Ecuaciones Lineales. Metodos Directos e Iterativosmetodos numericos itescoRombergmullerANALISIS+NUMERICOEjercicios Resueltos Interpolacion PolinomialAtenciónMatrizMÉTODOS NUMÉRICOS EN MATLABmetodos numericos para ingenieriaMÉTODO DE HORNERJacobi Program in FortranMétodos numéricos básicos para ingenieríahalleyeuler y heunMetodos de Solucion de Ecuaciones No LinealesMatrices de Diagonal Estrictamente DominanteMétodos NuméricosMétodo de SteffensenEjercicios Resueltos Integracion NumericaMÉTODOS DIRECTOS PARA SOLUCIÓN DE SISTEMAS ECUACIONES LINEALES.Una Introducción al Análisis NuméricoEjercicios resueltos Topología [E. Nieto]Newton RaphsonProblemas Resueltos de Analisis Numericos IMetodos Numericos Chapra 5 Edicionmatlab con metodos numericosMétodo de Newton-RaphsonMore from Daniel Nuñez HernandezRapid PCR Plastic Versus Glasslm352014_19_mercurio_5_1NCh0015-5-2000.pdfgcm-03.pdfFlujo de Fluidos en valvulas, - CRANE Co._1439.pdfNuevo documento de texto (2).txtOpcion 1Opcion 2Pieza TensorFlexion Wrist for MC Hand or ETDPHS 100-1A21ARTICULADOS- INTRODUCCION 2014Angular Contact Ball Bearings, Single RowDialnet IngenieriaInversaDeUnReductorDeTornilloSinfincoron 4271876 (1)Maxon Motor 118392ReductorMSDS-HPC Isometrica TuberiasTD11-03 Grafico de espacios maximos entre guias - (Macoga).pdfClave cengel tranfe.docx68917550-Guia-de-Valvulas.pdfFlujo compresible.pdf26_evaluacion.pdfAqui estaba la papa.pdfManualCompostajeparaMunicipios.pdfTest A la caldera Una vez finalizado el calorimetro.pdfManualCompostajeparaMunicipios.pdfManual-de-Compostaje-Casero-MMA.pdfỘģỘģ庨ħd - Analisis analitico.pdf

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