Source: https://www.scribd.com/document/52039575/Apuntes-de-analisis-numerico
Timestamp: 2016-08-27 05:17:20+00:00

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

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