Source: https://es.scribd.com/doc/52039575/Apuntes-de-analisis-numerico
Timestamp: 2016-02-06 18:49:45+00:00

Document:
SubirSign inJoinBooksAudiobooksComicsSheet MusicScribd Selects BooksHand-picked favorites from our editorsScribd Selects AudiobooksHand-picked favorites from our editorsScribd Selects ComicsHand-picked favorites from our editorsScribd Selects Sheet MusicHand-picked favorites from our editorsTop BooksWhat's trending, bestsellers, award-winners & moreTop AudiobooksWhat's trending, bestsellers, award-winners & moreTop ComicsWhat's trending, bestsellers, award-winners & moreTop Sheet MusicWhat's trending, bestsellers, award-winners & moreCategoriesArts & IdeasBiography & MemoirBusiness & LeadershipChildren'sComputers & TechnologyCooking & FoodCrafts & HobbiesFantasyFiction & LiteratureHappiness & Self-HelpHealth & WellnessHistoryHome & GardenHumorLGBTMystery, Thriller & CrimePolitics & EconomyReferenceReligionRomanceScience & NatureScience FictionSociety & CultureSports & AdventureTravelYoung AdultCategoriesArts & IdeasBiography & MemoirBusiness & LeadershipChildren'sComputers & TechnologyCooking & FoodFantasyFiction & LiteratureHappiness & Self-HelpHealth & WellnessHistoryHome & GardenHumorLGBTMystery, Thriller & CrimePolitics & EconomyReferenceReligionRomanceScience & NatureScience FictionSociety & CultureSports & AdventureTravelYoung AdultCategoriesAdaptationsChildren’sCrime & MysteryFictionHumorMangaNonfictionRomanceSciFi, Fantasy & HorrorSuperheroesYoung AdultPublishersArcanaArchie ComicsBOOM! StudiosDynamiteIDW PublishingKingstone ComicsMarvel ComicsSpace Goat ProductionsTop Cow ComicsTop Shelf ProductionsValiant Comics ZenescopeDifficultyBeginnerIntermediateAdvancedMixedInstrumentBrassDrums & PercussionGuitar, Bass, and FrettedPianoStringsVocalWoodwindsGenreClassicalCountryFolkJazz & BluesMovies & MusicalsPop & RockReligious & HolidayStandardsP. 1Apuntes+de+analisis+numericoApuntes+de+analisis+numerico|Views: 10.509|Likes: 151Publicado porDaniel Nuñez HernandezMore info:Published by: Daniel Nuñez Hernandez on Apr 01, 2011Copyright:Attribution Non-commercialAvailability:Read on Scribd mobile: iPhone, iPad and Android.download as PDF, TXT or read online from ScribdFlag for inappropriate content|Agregar a la colecciónSee moreSee lesshttps://es.scribd.com/doc/52039575/Apuntes-de-analisis-numerico06/04/2013pdftextoriginalSectionsINTRODUCCIÓNAritméticas de precisión ﬁnitaPráctica 1 (Aritméticas ﬁnitas, 2 horas)Fuentes de errores numéricosMétodo de la bisecciónMétodo de la Regula-falsi (regla de lo falso)Método de Newton-RaphsonEl método de la SecanteMétodo de MüllerPráctica 2 (Método de Müller, 4 horas)Cálculo de las raíces de un polinomioINTERPOLACIÓN DE FUNCIONES IInterpolación por polinomios de LagrangeImplementación de funciones elementalesANÁLISIS NUMÉRICO MATRICIAL IMétodo de CholeskyPráctica 4 (Método de Cholesky, 6 horas)Método de Crout para matrices tridiagonalesDiferenciación NuméricaDiscretización del LaplacianoDiscretización del gradienteIntegración NuméricaMétodos de Cuadratura de GaussANÁLISIS NUMÉRICO MATRICIAL IINormas de vectores y matricesCondicionamiento de una matrizCálculo de autovalores y autovectoresMétodo de la potencia inversaMétodo de Gauss-SeidelConvergencia de los métodos iterativosPráctica 7 (Método de relajación, 2 horas)INTERPOLACIÓN DE FUNCIONES IIInterpolación de HermiteInterpolación por splines cúbicosAproximación por mínimos cuadradosBIBLIOGRAFÍA BÁSICAANÁ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
Contenidos 1 INTRODUCCIÓN 2
2 ARITMÉTICAS DE PRECISIÓN FINITA Y FUENTES DE ERRORES NUMÉRICOS 2 2.1 Aritméticas de precisión ﬁnita . . . . . . . . 2 2.2 Práctica 1 (Aritméticas ﬁnitas, 2 horas) . 5 2.3 Fuentes de errores numéricos . . . . . . . . 7 3 CÁLCULO DE LOS CEROS DE UNA FUNCIÓN 8 3.1 Método de la bisección . . . . . . . . . . . . 8 3.2 Método de la Regula-falsi (regla de lo falso) 8 3.3 Método de Newton-Raphson . . . . . . . . . 8 3.4 El método de la Secante . . . . . . . . . . . 8 3.5 Método de Müller . . . . . . . . . . . . . . . 9 3.6 Práctica 2 (Método de Müller, 4 horas) . . 9 3.7 Cálculo de las raíces de un polinomio . . . . 10 3.7.1 Algoritmo de Horner para evaluar un polinomio en un punto . . . . . . 10 4 INTERPOLACIÓN DE FUNCIONES I 14 4.1 Interpolación por polinomios de Lagrange . 14 4.2 Error de interpolación de Lagrange y polinomios de Chebychev . . . . . . . . . . . . 15 4.3 Método de diferencias de Newton para el cálculo del polinomio interpolador de Lagrange . . . . . . . . . . . . . . . . . . . . . 15 4.4 Implementación de funciones elementales . . 18 4.4.1 Aproximación de la exponencial ex . 18 4.5 Práctica 3 (Aproximación de ex , 2 horas) . 18 4.5.1 Aproximación de funciones trigonométricas . . . . . . . . . . . . 18 4.5.2 Aproximación de la función ln(x) . . 19 5 ANÁLISIS NUMÉRICO MATRICIAL I 5.1 Método de Gauss . . . . . . . . . . . . . . . 5.2 Estimación del error de un método para resolver sistemas . . . . . . . . . . . . . . . . 5.3 Método de Cholesky . . . . . . . . . . . . . 5.4 Práctica 4 (Método de Cholesky, 6 horas) 5.5 Método de Crout para matrices tridiagonales 5.6 Subrutinas en Fortran 77 para la lectura y escritura en disco de vectores y matrices . . 19 19 21 21 22 22 23
6 DIFERENCIACIÓN E INTEGRACIÓN NUMÉRICA 6.1 Diferenciación Numérica . . . . . . . . . . . 6.2 Diferenciación numérica en dimensiones superiores . . . . . . . . . . . . . . . . . . . . 6.2.1 Discretización del Laplaciano . . . . 6.2.2 Discretización del gradiente . . . . . 6.3 Integración Numérica . . . . . . . . . . . . . 6.3.1 Métodos de Cuadratura de Gauss . . 6.3.2 Fórmulas de Integración Numérica Compuestas . . . . . . . . . . . . . . 6.4 Práctica 5 (Implementación Método de Integración de Simpson, 2 horas) . . . . . . . 6.5 Integración numérica en dimensiones superiores . . . . . . . . . . . . . . . . . . . . . . 7 ANÁLISIS NUMÉRICO MATRICIAL II 7.1 Normas de vectores y matrices . . . . . . . 7.2 Condicionamiento de una matriz . . . . . . 7.3 Cálculo de autovalores y autovectores . . . 7.3.1 Método de Jacobi . . . . . . . . . . 7.4 Práctica 6 (Método de Jacobi para el cálculo de autovalores y autovectores 6 horas) 7.4.1 Método de la potencia . . . . . . . . 7.4.2 Método de la potencia inversa . . . . 7.5 Métodos iterativos de resolución de sistemas lineales . . . . . . . . . . . . . . . . . 7.5.1 Método de Jacobi . . . . . . . . . . 7.5.2 Método de Gauss-Seidel . . . . . . . 7.5.3 Método de relajación . . . . . . . . . 7.5.4 Convergencia de los métodos iterativos 7.6 Práctica 7 (Método de relajación, 2 horas) 7.7 Método de Newton-Raphson para sistemas de ecuaciones no lineales . . . . . . . . . . .
24 24 25 26 26 27 27 28 29 29 31 31 33 33 34 36 36 37 38 39 39 40 41 42 42
8 INTERPOLACIÓN DE FUNCIONES II 43 8.1 Interpolación de Hermite . . . . . . . . . . . 43 8.2 Interpolación por splines cúbicos . . . . . . 43 8.3 La interpolación a través de la función seno cardinal . . . . . . . . . . . . . . . . . . . . 46 8.4 La interpolación a través de polinomios trigonométricos . . . . . . . . . . . . . . . . 46 8.5 Aproximación por mínimos cuadrados . . . 47 9 BIBLIOGRAFÍA BÁSICA 48
10 APÉNDICE A: Resumen de los comandos de UNIX 49 11 APÉNDICE B: Resumen del procesador de texto vi 49 12 APÉNDICE C: Algunos fallos comunes en Fortran 49 INTRODUCCIÓN El presente documento es un texto de referencia básico sobre los contenidos de la disciplina de Análisis Numérico en el contexto curricular de una Ingeniería Informática. Aunque el texto cubre los contenidos mínimos necesarios, resultará de gran interés para los alumnos complementar la información aquí suministrada con los textos de referencia básicos mencionados en la bibliografía. Muchas de las demostraciones de los resultados presentados se encuentran en este texto. En los casos en que las demostraciones no se incluyen, se suministra el libro y la página donde se encuentra tal demostración, para que el alumno interesado pueda estudiarla por su cuenta. En general, todos los temas presentados aparecen bien desarrollados en los libros de texto clásicos mencionados en la bibliografía. La única excepción es el tema de aritméticas de precisión ﬁnita, que se ha desarrollado en este texto con algo más de detalle y con un enfoque algo más moderno que en los libros clásicos, por considerar que, en el contexto de una Ingeniería Informática, este tema es de especial relevancia. El lenguaje de programación que se utilizará es el Fortran. Se ha elegido este lenguaje por ser la plataforma donde se han desarrollado habitualmente los grandes programas de cálculo numérico y por estar especialmente orientado al cálculo cientíﬁco. En el texto se va introduciendo este lenguaje de programación a través de programas ejemplo. Estos programas ejemplo se encuentran a disposición de los alumnos en el directorio de la asignatura /users/asignaturas/ii-an de la máquina serdis.dis.ulpgc.es. También se encuentra a disposición de los alumnos el ﬁchero an.h, donde se encuentran todas las subrutinas deﬁnidas en estos programas ejemplo. En el texto se proponen unas prácticas de laboratorio para realizar a lo largo de la asignatura. Para establecer el orden de impartición de los contenidos presentes en este documento se ha utilizado, como criterio preferente, la coordinación entre el programa de prácticas y el programa teórico de la asignatura, de tal forma que, con un desarrollo normal de la docencia, los contenidos teóricos sean presentados con antelación al desarrollo de las prácticas, comenzando las prácticas de laboratorio a partir de la segunda semana de clase. Para el buen seguimiento de la asignatura, resulta de gran interés tener cierta soltura en el manejo de los conceptos elementales del Análisis Matemático, el Álgebra, y la programación de Algoritmos. La materia expuesta en esta documentación está programada para ser impartida en un cuatrimestre a razón de 3 horas/semana en el aula y 2 horas/semana en el laboratorio informático, lo
que hace un total de, aproximadamente, 45 horas en aula (3 créditos teóricos) y 30 horas de laboratorio (2 créditos prácticos). Dado el escaso tiempo disponible, se han eliminado algunos temas clásicos de un curso completo anual de Análisis Numérico como son las ecuaciones diferenciales ordinarias y las ecuaciones en derivadas parciales. Normalmente, dichos temas se verán en detalle en asignaturas posteriores. Además, en lugar de presentar de forma exhaustiva todos los métodos numéricos que se pueden encontrar en los libros de Análisis Numérico clásicos, se ha optado por reducir los contenidos e impartir una selección de los métodos numéricos más representativos. ARITMÉTICAS DE PRECISIÓN FINITA Y FUENTES DE ERRORES NUMÉRICOS Aritméticas de precisión ﬁnita Un número entero z se representa en el ordenador a través de un número ﬁjo de bits (16 bits habitualmente), donde uno de los bits se utiliza para determinar el signo y los restantes para expresar el valor absoluto del número, de tal manera que la secuencia de bits a1 a2 a3 ......an donde ai = 0 o ai = 1, representa el valor absoluto del número | z |= an + an−1 2 + an−2 22 + ... + a1 2n−1 Así, utilizando 16 bits, el mayor número entero que podemos representar es 1 + 2 + 22 + ..... + 214 = 215 − 1 = 32767 Es decir, los número enteros que podemos expresar con una aritmética de 16 bits van desde −32767 hasta 32767. Para representar un número real y en el ordenador nos basaremos en el siguiente resultado: Teorema 1 Un número real positivo y se puede expresar como ∞ X an y = 2e 2n n=1
donde e es un número entero, a1 = 1, y para n > 1, an = 0 o an = 1. Demostración. Dado un número real positivo y, existe un entero e tal que 2e−1 ≤ y < 2e , y por tanto 2−1 ≤ y2−e < 1. Por otro lado, si deﬁnimos las sucesiones Sn y an de la siguiente forma: S1 = 2−1 , an = 1 y para n > 1 an = 0 si Sn−1 + an = 1 si Sn−1 + Sn =
1 2n 1 2n
> y2−e ≤ y2−e
entonces es claro que | Sn − y2−e |≤ 21 , y por tanto Sn → n y2−e lo que concluye la demostración del teorema.
n X ak
Ejemplo 1 Consideremos y = 10. 125, podemos expresar este número como 1 1 1 10. 125 = 24 ( + 3 + 7 ) 2 2 2 Es decir, e = 4, a1 = a3 = a7 = 1, y el resto de los an es 0. En este caso, el número de elementos an distintos de 0 es un número ﬁnito, en general no ocurre así. Evidentemente, cualquier número que tenga un número ﬁnito de elementos an distintos de 0 es un número racional y, por tanto, los números irracionales se representarán siempre con un número inﬁnito de elementos an no nulos. Sin embargo, como muestra el siguiente problema, existen otros muchos números además de los irracionales, que no se pueden representar con un número ﬁnito de elementos an no nulos. Problema 1 (2 puntos) Demostrar que al representar el número real 0.1 como 0.1 = 2e
∞ X an 2n n=1
Problema 4 (2 puntos) Calcular todos los números reales que se pueden construir tomando 5 bits de la forma siguiente: 1 bit para el signo, 2 bits para la mantisa (es decir t = 3, puesto que a1 = 1 y sólo se almacenan a2 y a3 ) y 2 bits para el exponente e, tomando como rango de e = −1, 0, 1, 2. Representar dichos números sobre una recta. Es importante resaltar que los números reales en una aritmética de precisión ﬁnita no están equiespaciados, es decir, los números están más cercanos entre sí cerca de 0, y más alejados al alejarnos de 0. En 1985, la sociedad I.E.E.E. presentó una serie de especiﬁcaciones estándares para la deﬁnición de una aritmética de precisión ﬁnita para los números reales. En este trabajo, se codiﬁca un número real en simple precisión utilizando 32 bits de memoria, de los cuales 23 bits se utilizan para la mantisa (es decir t = 24 puesto que a1 = 1 no se almacena), 1 bit se utiliza para el signo y 8 bits se utilizan para el exponente e, lo cual da un rango de 28 = 256 valores posibles para el exponente e. En este caso, se toma emin = −125 y emax = 128. Como puede observarse, el número total de exponentes posibles es 254, dos menos que los 256 posibles, ello se hace así, porque se reservan dos casos para tratar las denominadas excepciones, como se verá más adelante. Por tanto, el valor máximo que puede tomar un número real en esta aritmética es ymax = 2128 e
24 X 1 = 3. 4 × 1038 2n n=1
el número de elementos no nulos an es inﬁnito. Problema 2 (2 puntos) Representar el número 0.0 703 125 como 0.0 703 125 = 2e
Para deﬁnir una aritmética de precisión ﬁnita de número reales, lo que se hace habitualmente es discretizar la fórmula de representación anterior, tomando un número ﬁnito de valores posibles ai y un número ﬁnito de valores para el exponente e. Como puede observarse, cada valor ai viene representado por un bit. Además, puesto que el valor a1 es siempre igual a 1, no es necesario almacenar su valor en memoria al guardar un número real. Por tanto, en una aritmética de precisión ﬁnita, los números reales distintos de cero se representan como y = ±2e e
t X an 2n n=1
y el valor mínimo positivo es ymin = 2−125 e 1 = 1. 18 × 10−38 2
donde e varía entre dos valores limites emin ≤ e ≤ emax . Al valor t se le llama precisión de la aritmética. A la secuencia a1 a2 a3 ......at , (donde ai ∈ {0, 1}) se le denomina mantisa. Hay que hacer notar aquí que, dado que hemos impuesto siempre que a1 = 1, el número 0 debemos añadirlo a la aritmética, ya que 0 no se puede representar de la forma anterior. Problema 3 (1 punto) Calcular los valores positivos mínimo y máximo que puede tomar un número real en una aritmética de precisión ﬁnita en función de t, emin y emax .
Además, el número de combinaciones posibles que puede tener la mantisa es 224 ≈ 1. 68 × 107 . Es decir, la aritmética tiene una precisión de 7 dígitos decimales. Esta representación equivale a normalizar el número binario colocando la coma detrás del primer 1, representar en la mantisa la parte fraccionaria, puesto que la parte entera es un 1 implícito, y el exponente en exceso 127, es decir, sumando esta cantidad al exponente resultante de la normalización, de forma que el rango que va desde 1 hasta 254 representa los exponentes que van desde −126 hasta 127 (se reservan los valores 0 y 255 para las excepciones). También se deﬁne en este trabajo de I.E.E.E. un estándar para una aritmética en doble precisión. En este caso, se utilizan 64 bits para almacenar un número real, de los cuales 52 bits se utilizan para la mantisa (t = 53), 1 bit para el signo y 11 bits para el exponente, lo que da lugar a 211 = 2048 posibilidades de elección de exponente e. En este caso, se toma emin = −1021 y emax = 1024. Por tanto, el valor máximo que puede tomar un número real en esta aritmética es ymax = 21024 e
53 X 1 = 1. 78 × 10308 2n n=1
y el valor mínimo positivo es ymin = 2−1021 e 1 = 2. 23 × 10−308 2
Además, el número de combinaciones posibles que puede tener la mantisa es 253 ≈ 9. 0 × 1015 . Es decir, la aritmética tiene una precisión de 15 dígitos decimales. Evidentemente, estos estándares no se siguen al pie de la letra por los diferentes fabricantes de ordenadores. Así, aunque existe bastante homogeneidad en este sentido, los valores pueden cambiar ligeramente de una máquina a otra. Por otro lado, aunque en la mayoría de los ordenadores actuales la base de representación de los números es 2, todavía pueden encontrarse algunos sistemas donde la base es 10 ó 16. Nosotros no entraremos aquí a estudiar este tipo de bases. El estudio es básicamente el mismo, adaptándolo a la base de representación. Tratamiento de las excepciones en el estándar de I.E.E.E. Denominaremos excepciones a las expresiones que no se pueden expresar en una aritmética usual, como son el 0, √ el inﬁnito, operaciones no válidas (como −1), etc. Estas excepciones son tratadas en el estándar de I.E.E.E. de la siguiente forma: dentro de las posiciones de memoria dedicadas al exponente e de un número, se reservan dos, que corresponden a emin − 1 y emax + 1, para trabajar con con las excepciones. La regla que se utiliza es la siguiente: 1. Si el valor de una variable y tiene por exponente emax + 1 y todos los coeﬁcientes de la mantisa valen 0, entonces y se considera inﬁnito. Por ejemplo 1/0 debe dar inﬁnito. 2. Si el valor de una variable y tiene por exponente emax + 1 y algún coeﬁciente de la mantisa es distinto de 0, entonces y se considera que no √ un número es (NaN (Not a Number)). Por ejemplo −1 debe dar NaN. 3. Si el valor de una variable y tiene por exponente emin − 1 y todos los coeﬁcientes de la mantisa valen 0, entonces y se considera igual a 0. 4. Si el valor de una variable y tiene por exponente emin − 1 y algún coeﬁciente de la mantisa es distinto de 0, y se considera que no está normalizado (es decir a1 = 0) y el valor de y sería y=2
emin −1 t X an 2n n=2
A=A/2. IF(A.GT.0) THEN M=M+1 GOTO 1 ENDIF PRINT *,M END
Nota En Fortran 77 no es necesario declarar las variables. Por defecto, las variables cuyo nombre empieza por las letras I, J, K, L, M, N son variables enteras, y el resto son variables reales. Las cinco primeras columnas de cada línea de un programa Fortran 77 están reservadas para escribir un número de etiqueta. La columna 6 está reservada para indicar si una línea es continuación de la anterior (en este caso basta con escribir un carácter en la columna 6.) Un ∗ o una C en la primera columna indica que la línea es de comentario. En el Fortran 77 estándar no existe la instrucción W HILE. Sin embargo, como se muestra en el programa anterior, se puede simular fácilmente un W HILE con un GOT O. Los operadores de comparación en fortran 77 son: .GT., .GE., .EQ., .N E., .LT. y .LE., que signiﬁcan mayor que, mayor o igual a, igual a, no igual a, menor que y menor o igual a, respectivamente. Los operadores lógicos son .AN D. y .OR., que realizan las operaciones de conjunción y disyunción, respectivamente. Programa 2 Programa en Fortran 77 para calcular el mayor número positivo de una aritmética. El programa devuelve un entero M tal que 2M es el número cuya representación corresponde a la excepción que codiﬁca el inﬁnito. A=1. B=1. M=0 B=2.*A IF(B.GT.A) THEN A=B M=M+1 GOTO 1 ENDIF PRINT *,M END
Programa 1 Programa en Fortran 77 para calcular el menor número positivo de una aritmética. El programa devuelve un entero M tal que 2−M es el menor número real positivo normalizado. A=1. M=0
Problema 5 (2 puntos) Dada una aritmética de precisión ﬁnita cualquiera, calcular la distancia que hay entre el número 1 y su inmediato superior, es decir, el número que va después de 1, y la distancia entre el número 1 y su inmediato inferior. Vamos a llamar A, al conjunto de valores reales a los que da lugar una aritmética de precisión ﬁnita, es decir ( ) t X an e A = ±2 ∪ {0} 2n n=1 Dado un número real cualquiera y, al representarlo en una aritmética de precisión ﬁnita se produce un error de redondeo, llamaremos y ∈ A al número real que mejor e aproxima a y dentro de A. Deﬁnición 1 Dada una aritmética de precisión ﬁnita, se deﬁne la unidad de redondeo u como u = 2−t Por ejemplo, si t = 24 (reales en simple precisión) u = 2−24 = 5. 97 × 10−8 , y en doble precisión (t = 53), u = 2−53 = 1. 1 × 10−16 . Programa 3 Programa en Fortran 77 para calcular la unidad de redondeo de una aritmética. El programa devuelve un entero M tal que u = 2−M A=1. M=1 A=A*2. IF((1.+1./A).GT.1) THEN M=M+1 GOTO 1 ENDIF PRINT *,M END
serdis.dis.ulpgc.es y trabajar directamente sobre el terminal de conexión. Los ﬁcheros se pueden editar y corregir en entorno Windows. Compilar y ejecutar los programas 1, 2 y 3 para comprobar cuáles son el menor y el mayor número positivo, y la unidad de redondeo del ordenador en precisión simple. Dichos programas se encuentran en el directorio de la asignatura /users/asignaturas/ii-an de la máquina serdis.dis.ulpgc.es. Si ponemos en la cabecera del programa la instrucción IM P LICIT DOU BLE P RECISION (D)
cualquier variable cuyo nombre empiece por D será un número real en doble precisión. Hacer las modiﬁcaciones pertinentes en los programas 1, 2 y 3, para comprobar cuáles son el menor y el mayor número positivo, y la unidad de redondeo del ordenador en doble precisión. Además, √ hacer operaciones del tipo 1/0, 1/∞, ∞/0, ∞/∞, −1 e imprimir los resultados para ver cómo trata las excepciones el FORTRAN en la arquitectura de los ordenadores del laboratorio.
Problema 6 (4 puntos) Se considera una aritmética de 16 bits donde se dedica 1 bit al signo, 9 bits a la mantisa (t = 10) y 6 bits al exponente (emin = −30 emax = 31). Escribir, si es posible, los siguientes números en esta aritmética: 1. 2, y los números más cercanos a 2 por arriba y por debajo. 2. El cero, el inﬁnito y NaN. 3. Los números positivos mayor y menor de la aritmética (teniendo en cuenta las excepciones). 4.
Práctica 1 (Aritméticas ﬁnitas, 2 horas) La línea de comando para la compilación de un programa en Fortran 77 tiene la forma: > f 77 prog1.f − o prog1 Esta línea compila el programa prog1.f y genera el ejecutable prog1. Para compilar un programa en Fortran, necesitamos una máquina que tenga instalado el compilador. Por ejemplo, la máquina serdis.dis.ulpgc.es tiene dicho compilador. Para utilizar el compilador desde un entorno Windows, basta con conectarse a través de la utilidad SSH a
Problema 7 ¢ (2 ¡ 1 1 2 1 + 23 + 25 B 2 B+A y B−A
. puntos) ¡ Sean 1 = 23 1 + 26 + 2
= ¢A . Calcular
Problema 8 (2 puntos) Sean emin , emax , los valores mínimo y máximo del exponente e. Demostrar que si emin < e < emax , entonces los números: Ã t ! X an 1 ± t 2e 2n 2 n=1 pertenecen al conjunto A de números reales generados por la aritmética de precisión ﬁnita.
Un número real cualquiera z.LE. También se puede utilizar un criterio más simple.B.(TOL*(ABS(B)+10. salvo que A también sea 0.A continuación. 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.GT. Por ejemplo. Si un número real z veriﬁca que ymin <| z |< ymax .TOL) IF(ABS(A).**(10.**(10. como | A − B |≤| A | T OL pero en este caso le estamos dando una signiﬁcación especial a A con respecto a B. por tanto | z − z |≤| z | 2−t =| z | u e con lo que queda demostrado el teorema. los criterios anteriores quedan | A |≤| A | T OL lo cual es imposible (si T OL < 1).B. el test de parada incluye el hecho de que dos variables estén próximas entre sí.’A=B segun la tolerancia TOL’ STOP ELSE PRINT *. z2 ∈ A son distintos entonces e e Demostración: Ejercicio En muchos algoritmos. | B |} T OL Este criterio es simétrico en el sentido de que trata de igual modo los números A y B. Problema 9 (2 puntos) Dado un número z = e P an 2e t n=1 2n . | B |} + ) T OL
donde a0 = 1 y. con = 10−10 READ *. como a0 = 1.B. entonces e e donde z es el número más cercano a z en la aritmética. 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. en general. el número que está a la derecha de la desigualdad también pertenece a la aritmética de precisión ﬁnita. si B = 0. que tomaremos positivo sin pérdida de generalidad. B). se tiene que 2e < 2 | z | y. Para evitar este comportamiento.TOL). Teorema 2 Sean ymin . Además.ABS(B)) THEN IF(ABS(A-B).(TOL*(ABS(A)+10.0) THEN PRINT *.A.TOL IF(IGUAL(A. B son iguales con una tolerancia T OL (tomando el máximo de A.))) THEN IGUAL=0 RETURN ELSE IGUAL=1
Ahora bien. 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 |. Calcular el número inmediatamente inferior a él en dicha aritmética.’A distinto de B segun la tolerancia TOL’ STOP ENDIF END FUNCTION IGUAL(A.LE.EQ. | z2 |} u e e e e
. ymax los valores positivos menor e e y mayor de una aritmética de precisión ﬁnita. en una aritmética de precisión ﬁnita. Un resultado importante para la comparación de dos números es el siguiente: Teorema 3 Si z1 . Sea u la unidad de redondeo de dicha aritmética. 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 |. y por tanto | z − z |≤ e 2e 2−t 2
Programa 4 Programa en Fortran 77 que determina si dos variables A. se puede añadir al criterio un valor > 0 de la siguiente forma: | A − B |≤ (max {| A |. se puede expresar como ∞ X an z = 2e 2n n=1
Por el problema anterior.))) THEN IGUAL=0 RETURN ELSE IGUAL=1 RETURN ENDIF ELSE IF(ABS(A-B). an = 0 o an = 1. 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
las aritméticas estándares de ordenador trabajan en base 2.(1-D)*(10**10) END
Además.1 el ordenador. que son la suma. la resta de variables que tengan una magnitud cercana. podemos extraer que. Por ello. Hay que tener en cuenta que. este programa permite identiﬁcar la base de la aritmética con la que trabaja el ordenador. el resultado sea exactamente 1. Errores por Cancelación. números tan naturales para nosotros como 0. Solamente queremos mencionar que. para minimizar el efecto de los redondeos en las operaciones. el resultado sí es exactamente 1. Este resultado se pone de maniﬁesto en el siguiente programa Fortran: Programa 5 Programa en Fortran 77 para comprobar la diferencia entre trabajar en base 10 y trabajar en base 2. Este tipo de errores se produce al realizar un cambio de base para representar un número real. al realizar operaciones sobre una variable. se deberá evitar. el resultado se redondea para pasarlo a la precisión inicial.2**7 B=B+A CONTINUE DO 2 K=1. Como vimos en la sección anterior. 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. a menudo. en la medida de lo posible. quedando la aportación de los dígitos de menos valor. si sumamos 128 = 27 veces el número 2−7 .01. para ser más precisos numéricamente. destacaremos 3 tipos: Errores de redondeo. este error está controlado por la denominada unidad de redondeo. ﬁnalmente. Estos errores se producen al restar números de aproximadamente la misma magnitud. En los algoritmos. al tomar un número real z y aproximarlo en la aritmética por el valor z ∈ A más próximo. Son los que se producen al ”redondear” un número real para poder expresarlo en una aritmética de precisión ﬁnita. Nosotros no vamos a entrar en este curso en cómo se pueden deﬁnir algorítmicamente estas operaciones. en el programa Fortran anterior. 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 . 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. Sin embargo.1 no pueden representarse de forma exacta en una aritmética en base 2. cuando sumamos 100 veces el número 0. al representar 0.
Fuentes de errores numéricos Dentro de las posibles fuentes de errores numéricos. Problema 10 (1 punto) Calcular las raíces del polinomio P (x) = x2 − 2x + 0. existen 4 operaciones básicas.) B=0
. los humanos pensamos y razonamos en términos de números en base 10. Por ejemplo. se cancelan las partes signiﬁcativas. los errores de redondeo se van acumulando en la parte menos signiﬁcativa del número (los dígitos de menos valor).01 evitando los errores de cancelación. que es donde más error hay. 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. se ha utilizado este fenómeno de cancelación para poner de maniﬁesto la diferencia entre trabajar con bases distintas. Por lo tanto. que corresponde a los dígitos de mayor valor.(1-B)*(10**10) PRINT *. realizar la operación en una aritmética de mayor precisión y. Esto quiere decir que. Por ejemplo. en los algoritmos. pero.100 D=D+C CONTINUE PRINT *. parece razonable pensar que. no es así. muchas veces se intenta evitar la posibilidad de restar 2 números que pudieran ser de magnitud parecida. al restar dos números de magnitud parecida. y este pequeño error de redondeo se puede ir propagando hasta producir errores apreciables.01 D=0 DO 1 K=1. A=2**(-7. a continuación. u = 2−t . en la conocida fórmula del cálculo de raíces de un polinomio de grado 2. va a producir un pequeño redondeo. de tal forma que. la multiplicación y la división de números reales dentro de la aritmética. la resta. Por ejemplo. que es como solemos hacerlo. Sin embargo.
C=0. Como conclusión de este apartado. el error de redondeo tiene la e expresión: | z − z |≤| z | u e Errores de cambio de base. Por ejemplo. Como vimos en la sección anterior.RETURN ENDIF ENDIF END 1 Asociado a cualquier aritmética de precisión ﬁnita de números reales. dejando relativamente intacta la parte más signiﬁcativa del número.
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.TOL.0) THEN PRINT *. los valores de x para los cuales f (x) = 0. 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. 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. entonces f (x) = x2 − A = 0. es decir f (x) ≈ f (x0 ) + f 0 (x0 )(x − x0 ) que corresponde a un polinomio de grado 1. sin duda. Dada una aproximación inicial de la raíz x0 . a partir de x0 . f (a) · f ( a+b ) < 0.’El numero A no es positivo’ STOP ENDIF X0=(1+A)/2.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. Problema 15 (3 puntos) Escribir el código en Fortran 77 para implementar el método de la Regula-falsi. uno de los métodos más importantes y útiles para el cálculo de raíces. READ *. b]. y a continuación se calcula x1 como el cero de este polinomio.A. b] por la mitad de la siguiente forma: Se toma el punto medio a+b . 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. 2 si f ( a+b ) · f (b) < 0 entonces hacemos a = a+b y volvemos 2 2 a subdividir el nuevo intervalo [a. de forma general. esto es. En caso contrario. por el contrario.*X0) IF(IGUAL(X0.
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. a partir de x0 una secuencia xn de valores que van aproximando la raíz. 0]. b] donde la función f (x) cambia de signo.
Método de la bisección Se considera un intervalo [a. es decir: f (x0 ) x1 = x0 − 0 f (x0 ) y por tanto. entonces hacemos b = a+b y volvemos 2 2 a empezar. El método consiste en ir dividiendo el intervalo [a. f (x).X1. Es decir. es decir f (a)·f (b) < 0. en el razonamiento anterior. de teniendo en cuenta que si x = A.
El método de la Secante Este método es una variante del método de Newton para el caso en que no sea posible calcular la derivada de f (x)
. 2]. f (b)) con el eje x.
Método de Newton-Raphson Éste es.X0 STOP ELSE X0=X1 ENDIF CONTINUE PRINT *. obtenemos.TOL). 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.’LA RAIZ DE A ES’. y un número máximo de iteraciones N max .Nmax IF(A. Si. Si f ( a+b ) = 0 ya 2 2 hemos encontrado la raíz x = a+b . b] van aproximando la raíz. se considera el punto de intersección de la recta que pasa por los puntos (a. una aproximación mejor x1 de la raíz. Las sucesivas subdivisiones del intervalo [a.Nmax X1=X0-(X0*X0-A)/(2.’No máximo de iterac.EQ. se busca.LE. 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. f (a)) y (b. DO 1 K=1.0) THEN PRINT *.
de una forma analítica. calculamos xn por el método de Newton-Raphson. en el sentido de que. para determinar la igualdad entre dos números. 0). x2 = (1.1. x0 . 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. T OL. f (x) = x2 + 1 2. Utilizar el método para calcular los posibles ceros de las siguientes funciones: 1.
f 0 (xn−1 ) ≈
f (xn−1 )−f (xn−2 ) xn−1 −xn−2
f 00 (xn−1 ) (xn−1 2
− xn−2 )
Como veremos posteriormente. 0). x1 y x2 . el número máximo de iteraciones. 0). x1 = (2. Para el ejemplo 5 tomar como datos iniciales x0 = (3. 0). 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). Dicha raíz será la aproximación xn de la raíz de f (x) en la etapa n. x1 = (2. 0). se sustituye el valor f 0 (xn ) en el algoritmo. Para el ejemplo 3 tomar como datos iniciales x0 = (3. f (xn−2 )) y (xn−1 . La función a la que se le calculan los ceros se deﬁne en el propio cuerpo del programa.
xn−3 −xn−1
. x2 = (1. x1 = (2. 0). (xn−2 . x1 = 1. x1 = (2. 0). Crear un programa en Fortran 77 que tenga como datos de entrada: las tres primeras aproximaciones de la raíz. Para el ejemplo 1 tomar como datos iniciales x0 = (3. En el caso en que f 00 (xn−1 ) = 0. x1 = (2. f (xn−1 )). f (x) = (x2 + 1)x 3. x0 y x1 . x1 = (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. Para iniciar el algoritmo. y calcular posteriormente las derivadas de dicha parábola. 0).CY END FUNCTION CF(CX) IMPLICIT COMPLEX(C) CF=SQRT(CX) END
Práctica 2 (Método de Müller. Es una generalización del método de Newton-Raphson. IMPLICIT COMPLEX (C) CX=(-1. N max. 0). 0). 4 horas) Implementar el método de Müller. Programa 7 Programa en Fortran 77 donde se muestra un ejemplo de manejo de números complejos. 0). y la tolerancia T OL. la elección de las fórmulas anteriores equivale a aproximar f (x) por la parábola que pasa por los puntos (xn−3 . Problema 18 (3 puntos) Escribir un programa en Fortran 77 que implemente el método de la Secante utilizando reales de doble precisión. x2 = (1. nos quedamos con los términos hasta el orden 2. para determinar la igualdad de dos números. Problema 16 (1 punto) Calcular una iteración del método de Newton-Raphson para calcular un cero de la función f (x) = x3 − 3 partiendo de x0 = 1. 0). f (x) = x − 2 5. f (x) = ex − 1 4. son necesarias dos aproximaciones iniciales. En este caso. 0) y x0 = (1. 0).0) CY=CF(CX) PRINT *. el número máximo de iteraciones N max. Para obtener una aproximación xn mejor de la raíz calculamos los ceros del polinomio de segundo grado anterior. 0). Los datos de entrada son las aproximaciones iniciales. como por ejemplo polinomios. x2 = (1. x0 y x1. En el caso en que no conozcamos analíticamente el valor de la primera y segunda derivada de f (x). f (x) = 1
Método de Müller Este método es de utilidad para calcular raíces complejas de funciones. f (xn−3 )) . x2 = (0. 0). por el valor f (xn ) − f (xn−1 ) xn − xn−1 que corresponde a una aproximación de f 0 (xn ). podemos utilizar las siguientes aproximaciones: f 00 (xn−1 ) ≈ 2
f (xn−2 )−f (xn−3 ) f (xn−1 )−f (xn−2 ) − xn−2 −xn−3 xn−1 −xn−2
Nota: Utilizar como tolerancia T OL = 0. 0). nos quedamos con aquélla que sea más cercana a xn−1 . y la tolerancia.0001 y N max = 100. 0) x2 = (1. en lugar de quedarnos con la parte lineal del desarrollo de Taylor de la función.01. Para el ejemplo 2 tomar como datos iniciales x0 = (3. Para el ejemplo 4 tomar como datos iniciales x0 = (3.
Este teorema permite calcular el polinomio y su derivada en un punto de forma muy sencilla.. almacenándolos en las variables P X y P P X..487%.. 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. Este ejemplo muestra como un problema ﬁnanciero sencillo nos lleva a la necesidad de calcular los ceros de un polinomio.. + (b0 − b1 x0 ) Por último..’EscriBir Coef. aportación que se va incrementando cada año en un 10%. PARAMETER(NMAX=1000) DIMENSION A(0:NMAX) COMMON/POL/PX. Además. con frecuencia. Sea P (x) = an xn + an−1 xn−1 + .. Ejemplo 2 Actualmente están muy de moda los planes de pensiones. entonces. por su gran utilidad. A menudo. visto anteriormente.
Algoritmo de Horner para evaluar un polinomio en un punto Dado un polinomio P (x) = an xn + an−1 xn−1 + .. éste se puede expresar también de la forma siguiente: P (x) = a0 + x (a1 + x (a2 + x (a3 + x(. aparecen problemas que. para calcular i. por ejemplo. + b1 .Cálculo de las raíces de un polinomio Los polinomios son un tipo particular de funciones que. dado que ak = bk − bk+1 x0 y an = bn . el método de Müller.’ DO 1 K=0. + a0 ...PPX END
(100. Programa 8 El siguiente programa en Fortran 77 calcula la evaluación de un polinomio y su derivada en un punto X. 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 + .000 pesetas todos los años. aunque también hay que indicar que existen algoritmos versátiles para el cálculo de las raíces complejas..X) PRINT *. obtenemos la igualdad anterior teniendo en cuenta que (x − x0 )Q(x) + b0 = bn xn + (bn−1 − bn x0 )xn + ...PPX PRINT *.’Escribir Grado del Polinomio’ READ *..1) (1.N IF(N. ’Grado Superior al Maximo’ STOP ENDIF PRINT *.PX PRINT *.NMAX) THEN PRINT *.N READ *. por ejemplo: si usted aporta durante 30 años 100.GT. Mi experiencia como docente en esta disciplina es que. requieren el uso de alguna de las técnicas presentadas en esta asignatura.. para su resolución. + a0 . Ahora bien. + b1 0 0 Demostración Sea el polinomio Q(x) = bn xn−1 + bn−1 xn−2 + ..000. Si llamamos i al interés nominal anual que se aplica al dinero..’Escribir valor de X’ READ *.A. necesitamos evaluar tanto el polinomio como su derivada.000
Ahora bien. una vez terminada la carrera y en el desarrollo de la actividad profesional..000) (1. El siguiente
... debemos calcular las raíces del polinomio en i dado por P (i) =
(100.000. + x (an−1 + xan ))))) . 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. requieren un análisis algo más detallado. 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. es decir el primer año 100..000 de pesetas. Teorema 4 (Método de Horner). obtenemos P 0 (x) = (x − x0 )Q0 (x) + Q(x) de donde sale obviamente que P 0 (x0 ) = Q(x0 ).. + i)30−n = 26.000.000) (1..000.. etc.A(K) PRINT *.000
El cálculo de las raíces de este polinomio nos lleva a i = 4. + i)30−n − 26. Polin.’P ‘(X)= ’. Nos ocuparemos sólo de las raíces reales de los polinomios.’P(X)= ’. El siguiente ejemplo es una buena prueba de ello. Las entidades ﬁnancieras venden a sus clientes los planes de pensiones de la siguiente forma. si queremos utilizar un método de cálculo de raíces como el de Newton-Raphson..1) (1..000.. el segundo año 110..X CALL HORNER(N... le aseguramos que al ﬁnal del trigésimo año tendrá a su disposición la cantidad de 26. como. como muestra el siguiente programa Fortran. Veamos que se veriﬁca que P (x) = (x − x0 )Q(x) + b0 Efectivamente..
...... hay un único cambio de signo y hay una raíz positiva.-1 PX=PX*X+A(K) PPX=PPX*X+PX CONTINUE PX=PX*X+A(0) END
Demostración [Is-Ke] Pg.. Por tanto..X) DIMENSION A(0:*) COMMON/POL/PX... Efectivamente.. entonces las m raíces distintas |an | x1 < x2 < . debe incluir en su inicio la misma sentencia COM M ON. entonces el número de raíces positivas es igual al número de cambios de signo en los coeﬁcientes an .. En este caso. o bien ese mismo número menos un número par..A.. 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 + . de tamaño N M AX + 1..n−1 |ak |) = >0 |x| − 1 = |an | |x|n − max |ak |
n−1 X k=0
Demostración Es inmediato. ≤ x0 1 2 m−1 ≤ xm ≤ Pmax Volviendo a aplicar este razonamiento sucesivamente sobre P 0 (x). Teorema 8 La derivada k − esima P k) (x) del polinomio ´ P (x) = an xn + an−1 xn−1 + . Para la estimación del número de raíces reales negativas...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. y localizarlas en un intervalo.n−1 1 − |x| n |x| n = ≥ |an | |x| − max |ak | k=0...n−1 maxk=0... < xm de P (x) están intercaladas con las 0 raíces x0 < x0 < . se aplica el teorema anterior cambiando x por −x... P 00 (x). Ejemplo 3 Sea P (x) = 3x4 + 10x3 − 10x − 3. 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 + .n−1 | ak | −1 − .. es decir −Pmax ≤ x1 ≤ x0 ≤ x2 ≤ x0 ≤ .... 126. Los dos resultados anteriores permiten aislar las posibles raíces de P (x) de la forma siguiente: Si llamamos maxk=0. y numerados desde 0 hasta N M AX. Por tanto. para intercalar los ceros de una derivada con los ceros de la siguiente. Si cambiamos x por −x... La declaración COM M ON/P OL/P X...
Nota: La declaración P ARAM ET ER(N M AX = 1000) permite deﬁnir constantes. P P X..
.. + a0 con an 6= 0. los signos de los coeﬁcientes son: + + −−. − 1 . . |P (x)| ≥ |an xn | −
k=0. utilizando el algoritmo de Horner.. evaluar el polinomio y su derivada en el punto x = 2..1
SUBROUTINE HORNER(N.1. etc.n−1 | ak | maxk=0... Teorema 7 Entre dos raíces de una función derivable f (x) hay una raíz de f 0 (x). Para que una subrutina pueda hacer uso de esas variables..n−1 |x| − 1 |x|n (|an | (|x| − 1) − maxk=0... + a0 . −3..1 + | an | | an | Demostración Veamos que si |x| > 1 + entonces |P (x)| > 0. de reales en precisión simple. 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 .
Problema 21 (1 punto) Calcular el número máximo de raíces positivas y negativas del polinomio x5 −35x3 +30x2 + 124x − 120... La declaración DIM EN SION A(0 : N M AX) deﬁne un vector A.+ak x x (n − k)! (n − k − 1)! 1
1 − |x| ≥ k=0. < x0 1 2 m−1 de P (x). los signos de los coeﬁcientes son + − +−.n−1 |ak | . Demostración Teorema de Rolle.. a0 (saltando los posibles coeﬁcientes nulos)..n−1 |ak | Pmax a 1 + . las raíces son x = 1. entonces las raíces reales de P (x) están en el intervalo ∙ ¸ maxk=0.. derivando sucesivamente el polinomio P (x)... −1.a0 es P k) (x) = k! an n! n−k an−1 (n − 1) ! n−k−1 + +. P P X deﬁne la zona de memoria denomina P OL donde se encuentran las variables globales P X. hay 3 cambios de signo y hay una o tres raíces negativas..
5 x -2.166. Existen métodos mejores para el cálculo de raíces de polinomios. 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. .. Se calcula la raíz x1 del Polinomio P n−1) (x) (que es un polinomio de grado 1) 3. <
xk+1) mk+1
< Pmax
Al ﬁnal del procedimiento.25 2.25 -25 25 0 0 1. Para k = n − 2. 1. Pmax ]
2. y cuya gráﬁca es
20 50 0 -2. 1 −Pmax <
k+1) x1
-2. cuya raíz es x = 0.858 para el inter12
. 30! = 2. −2.25 2. 574. 7.25 0 1.5 -1.5 -1. Por ejemplo.858. utilizando cualquier método numérico de los vistos anterioremente.cuyas raíces son x = −0.25 0 1. 2.25
50 25 0 0 -25 -50 -75 -100 1. el método de la Regula-falsi.25 2.5 -1. pero que utilizan técnicas más complejas.. 3. 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).5 x
Se calculan las raíces de P k (x) en los intervalos <
k+1) x2
< . funciona razonablemente bien para grados de polinomios pequeños..1. El método presente en el siguiente programa.25. 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.cuyas raíces son x = −1. Ejemplo 4 Consideremos el polinomio P (x) = x4 − x3 − 7x2 +x+6. 6 × 1032 .5 -1. que tiene por raices x = 1. buscamos las raíces de P 00 (x) en esos intervalos. 358 y cuya gráﬁca es
25 0 -2.5 x
Polinomio P (x) = x4 − x3 − 7x2 + x + 6
-50 -75
La derivada de este polinomio es P 0 (x) = 4x3 − 3x2 − 14x + 1.166]. habremos aislado completamente a las raíces de P (x). es decir x = 0. 2. obteniendo −0. por ejemplo. Puesto que hay cambio de signo de P 00 (x) en cada uno de estos intervalos.25. Por tanto.25] y [0.5 x
Polinomio P 00 (x) = 12x2 − 6x − 14
La derivada tercera de este polinomio es P 000 (x) = 24x − 6. puesto que su utilización requiere el cálculo de factoriales. tenemos que Pmax = 8. 05 × 10−2 ..25 2. las raíces están en el intervalo [−8. 0. 8]. −1. Este procedimiento se puede utilizar para grados relativamente pequeños (n < 30).25. que se dispara rápidamente.. por tanto los ceros de P 00 (x) estarían en los intervalos [−2. Se parte del intervalo [−Pmax . 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. Para este polinomio.
PARAMETER(NMAX=30) DIMENSION A(0:NMAX). 253. tiene’. 7. Buscamos.5. [−0. N maxx).N. 7. 358 para el intervalo [0.858.N IF(R(K-1). Polin.N 7 PI(K)=PMAX DO 4 K=N-2.25] y 1.A. donde están los coeﬁcientes del polinomio. También se deﬁne la función auxiliar RP (N. 574]. 1.Nmaxx. 2.5].M-1 9 PRINT *. F(0:NMAX).R. tomando como valor inicial el punto medio del intervalo [X1.X2) THEN RP=X1 RETURN ENDIF RP=(X1+X2)/2. 253] y [2. 1.574. Dicha subrutina tiene como parámetros un vector A().166].R(K) END
DIMENSION A(0:*).EQ.TOL. 05×10−2 ]. X1.EQ. DO 2 K=1. N maxx.TOL.Nmaxx) PRINT *. para NewtonRaphson’ READ *.253.’Escribir Grado del Polinomio’ READ *.R.N-K PI(L)=RP(N-K. −0.PI(L).A. obteniendo x = −1. 3. 4.PPX R(L)=1.PPX **** Calculo de los factoriales F(0)=1.LT.L-1) 6 CONTINUE 4 CONTINUE *** Pasamos las raices al vector R() M=0 DO 8 K=1. ’Escribir Tolerancia’ READ *. [−1.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. 0.N 10 PI(2)=PMAX *** Calculo de los coeﬁcientes del *** polinomio derivada DO 7 K=2. 2. 05 × 10−2 y 2.X2.0.M. DO 1 K=1. T OL. Max. R(0:*). N. [7. Por tanto.N 1 READ *. A.0. ’Grado Superior al Maximo’ STOP ENDIF PRINT *. con la que consideramos que dos números son iguales.EQ. PI(0)=-PMAX PI(1)=-(A(N-1)*F(N-1))/(A(N)*F(N)) DO 10 K=2.Nmaxx) PARAMETER(NMAX=30)
FUNCTION RP(N.Nmaxx.R.A(K) PRINT *. AP(0:NMAX). un vector R().RP) IF (PPX.N IF(N. Problema 22 (2 puntos) Aislar en intervalos las raíces del polinomio P (x) = 20x3 − 45x2 + 30x − 1.ABS(A(K)) THEN PMAX=ABS(A(K) ENDIF 3 CONTINUE PMAX=PMAX/ABS(A(N))+1.EQ.R(0:NMAX-1) COMMON/POL/PX.TOL PRINT *.0) THEN R(M)=PI(K) M=M+1 ENDIF 8 CONTINUE ICEROPOL=M END
FUNCTION ICEROPOL(A.TOL. ﬁnalmente.358.Nmaxx M=ICEROPOL(A.358] y [1.’Escribir Coef. −1. que devuelve las raíces reales de un polinomio. que devuelve la raíz del polinomio que se obtiene aplicando el método de Newton-Raphson. para el proceso de Newton-Raphson.0.) THEN IF(PX.-1 DO 5 L=0.
Programa 9 Programa en Fortran 77 donde se implementa la función ICEROP OL(A. iter.25.TOL. X2]. R. −1.N.’El Pol.N-1 IF(PMAX. Buscamos ahora las raíces de P 0 (x) es esos intervalos.’ raices’ DO 9 K=0. 05×10−2 .R. X2. L). el grado del polinomio N.X1. 8]. donde se guardan las raíces del polinomio una vez calculadas. IF (X1.’ DO 1 K=0.L) DIMENSION A(0:*). la tolerancia T OL.) THEN
. los posibles ceros de P (x) estarán en los intervalos [−8.AP.NMAX) THEN PRINT *.N 2 F(K)=F(K-1)*K *** Calculo intervalo inicial PMAX=ABS(A(0)) DO 3 K=1.166. ’Escribir No.GT. Por tanto. las raíces de P (x) en cada un de esos intervalos y obtenemos x = −2. T OL.R(0:*) COMMON/POL/PX.Nmaxx CALL HORNER(N. R. PI(0:NMAX+1) COMMON/POL/PX.PI(L1).858].PPX PRINT *.valo [−2. las posibles raíces de P 0 (x) estarán en los intervalos [−4. y el número máximo de iteraciones N max xx. 574.
...5 -1 -0. 2 2
ΠN i (xi − xj ) j6=
estos polinomios base tienen la propiedad fundamental siguiente ½ 1 si i = j P i (xj ) = 0 si i 6= j Por tanto. es el único polinomio de grado menor o igual que N tal que PN (xi ) = f (xi ) ∀i = 0. calcularíamos los polinomios base: (x + 1)(x − 1) −1 x(x − 1) 2 x(x + 1) 2
P 0 (x) = P 1 (x) = P 2 (x) =
siendo el polinomio interpolador:
P2 (x) = e0 Problema 23 (2 puntos) Aislar en intervalos las raíces del polinomio P (x) = 2x3 + 3x2 − 12x + 1. π .1
R(L)=0.5
1.. N.N . el polinomio interpolador de Lagrange en estos puntos.RP. a partir del conocimiento del valor de una función (y eventualmente de sus derivadas) en un conjunto ﬁnito de puntos.
En la siguiente ﬁgura comparamos la gráﬁca del polinomio P2 (x) (trazo continuo) con la gráﬁca de la función ex (trazo discontinuo)
Interpolación por polinomios de Lagrange Sea una función f (x) que conocemos en un conjunto ﬁnito de valores {xi }i=0.. RETURN ELSE RETURN ENDIF ELSE RP1=RP-PX/PPX IF(IGUAL(RP1. Por tanto Q(x) ≡ 0 y P (x) = PN (x). por tanto. vamos a interpolarla en los puntos x0 = 0. 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. Entonces.
PN (x) se puede expresar en término de los denominados polinomios base de Lagrange P i (x). aproximar el valor de la función fuera de ese conjunto ﬁnito de puntos.. . El polinomio interpolador de Lagrange PN (x) de f (x) en los puntos {xi }i=0. sabemos que f (xi ) = fi . RETURN ELSE RP=RP1 ENDIF ENDIF CONTINUE END
Ejemplo 5 Consideremos una función f (x) = ex . Para calcular P2 (x). π y 3π .. Es decir.0) THEN RP=RP1 R(L)=0... posee N + 1 raíces.EQ. x1 = −1 y x2 = 1.TOL)..5
x(x + 1) (x + 1)(x − 1) x(x − 1) + e−1 +e −1 2 2
INTERPOLACIÓN DE FUNCIONES I El problema general de la interpolación de funciones consiste en. N
-1. salvo que Q(x) sea identicamente igual a cero. . 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. el polinomio interpolador de Lagrange puede expresarse como PN (x) =
Teorema 9 El polinomio interpolador de Lagrange es el único polinomio de grado igual o inferior a N tal que PN (xi ) = f (xi ) ∀i = 0... N
f (xi )P i (x)
Demostración Sea P (x) un polinomio de grado inferior o igual a N que veriﬁque que P (xi ) = f (xi ) ∀i = 0.N . lo cual es imposible.
b]. b] se obtiene fácilmente transformando el intervalo [−1.N ⊂ [a. entonces f (x) − PN (x) = f N +1) (ξ) N Π (x − xi ) (N + 1)! i=0
entonces max | ≤ ΠN (x i=0
x∈[a. que viene determinado por el siguiente teorema. Por tanto. 1] se encuentra en [Ki-Ch] Pg. b].. N xi = a + 2 2N + 2
Ejemplo 6 Se considera [a. Los puntos de interpolación dados por el teorema anterior son: x0 x1 x2 x3 x4 x5 = = = = = = . 982 96 . b]. π . b].. b] = [−1. 1]. 629 41 . Un método más directo para el cálculo de PN (x) es el denominado método de diferencias de
. 1] en [a. 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. su función derivada φ0 (t) tiene al menos n ceros repartidos entre los ceros de φ(t).b]
− xi ) |=
Demostración La demostración para el intervalo [−1. debemos cambiar todos los polinomios base de Lagrange. 1] al interpolar la función cos(x) en los puntos descritos en el ejemplo anterior. Para ello. b] Se consideran los puntos xi dados por µ µ ¶¶ b−a 2i + 1 1 + cos π i = 0. φ00 (t) tiene al menos n − 1 ceros y así sucesivamente hasta llegar a φN+1 (t). al aproximar f (x) por el polinomio interpolador PN (x) en un intervalo [a. que tiene al menos 1 cero. Si llamamos ξ a dicho cero. los valores óptimos de interpolación xi dados por la fórmula anterior son las raíces de los denominados polinomios de Chebychev. 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. y PN (x) su polinomio interpolador de Lagrange en los puntos {xi }i=0. utilizando este resultado. TN (x). elegiremos los puntos xi tales que ΠN (x − xi ) sea lo más pequeño posible en i=0 [a. el error de interpolación máximo viene determinado por: | f (x) − PN (x) |≤ maxx∈[a. el error de interpolación es cero y por tanto la fórmula anterior es válida.b]
x∈[a. 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. 2π] interpolando en los puntos 0. despejando y sustituyendo λ por su valor. cómo elegirlos de tal forma que el error de interpolación sea mínimo. 292-294. 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. π y 2 3π 2 .Error de interpolación de Lagrange y polinomios de Chebychev Evidentemente. Por tanto. y que nosotros podamos elegir los valores de interpolación xi . b]. Análogamente. 703 7 × 10−2
Problema 26 (2 puntos) Calcular el error máximo de interpolación en el intervalo [0. Demostración Si x = xi . y acotarlo superiormente. b] se comete. En el caso de que [a. Además. si queremos añadir un nuevo punto de interpolación. La cuestión que vamos a abordar en este apartado es. un error de interpolación. 370 59 ... 853 55 . .. e
max | ΠN (x − xj ) | e j=0
1 ≤ 2N
donde ξ es un valor intermedio perteneciente a [a. 146 45 1. y un intervalo [a. b] y x ∈ [a. b] = [0. 1] y N = 5 (es decir 6 puntos de interpolación). La demostración para un intervalo cualquiera [a. obtenemos φN+1 (ξ) = f N+1) (ξ) − λ(N + 1)! de donde. en el caso en que queramos interpolar una función en un intervalo [a. Teorema 10 Sea f (x) una función. en general.. Teorema 11 Sea N ≥ 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 . obtenemos el resultado del Teorema.
los coeﬁcientes que acompañan a la potencia xk−1 en ambos polinomios coinciden y... 2] = 2 e3 − 2e2 + e1 f [1.. .. ...... xi+k ] = xi+k − xi
...... el polinomio P2 (x) lo expresamos como P2 (x) = 1 + (e − 1)x + e2 − 2e + 1 x(x − 1) 2
xi+j = bk−1 − bk
Despejando obtenemos ak = bk−1 − ak−1 xk+i − xi
Finalmente obtenemos el resultado del teorema.. . xi+k ] indica. x3 = 3.. f [xi+1 . despejando obtenemos e2 − P1 (2) a2 = 2 Por tanto. veriﬁcan las siguientes propiedades: f [xi ] = f (xi ) f [xi+1 ] − f [xi ] f [xi . xi+k ]
Como veremos en el teorema siguiente. xi+k ] = f [xi+k . Como el polinomio interpolador es único. que se denominan diferencias divididas de Newton.. . es decir. xi+k y.. xi+k−1 ] f [xi ... xi+k .... 1] = e1 − 1 f [1. xi+k ] − f [xi ... x1 = 1 y x2 = 2. .. . x1 = 1. xi . ... . xi+k ] no depende del orden en que tomemos los puntos xi . xi ] Consideremos ahora el polinomio interpolador Qk (x) que interpola en los puntos xi+k ... . .. xi+1 ] = xi+1 − xi . xi+k−1 ] = f [xi+k . 1. xi+k−j ] Por la unicidad del polinomio interpolador obtenemos que Pk (x) = Qk (x) y..... 2] = e2 − e1 f [2. los coeﬁcientes f [x0 . xi+1 ] = f [xi+1 .. teniendo en cuenta que ak−1 bk−1 = f [xi .. despejando obtenemos a1 = e − 1 Por último P2 (x) = P1 (x) + a2 x(x − 1)
Demostración En primer lugar. f [xi .Newton.. 2. Qk (x) se puede escribir como Qk (x) = b0 + b1 (x − xi+k ) + b2 (x − xk+i )(x − xk+i−1 ) + .. El método consiste en ir calculando progresivamente los polinomios Pk (x) que interpolan la función en los puntos x0 . observamos que f [xi ... . xk ] Ejemplo 7 Vamos a interpolar la función f (x) = ex en los puntos x0 = 0. . obtenemos el polinomio interpolador de la siguiente forma: f [0. . el coeﬁciente que acompaña a la potencia xk en el polinomio interpolador Pk (x) para los puntos xi .. para cada Pk (x). donde bj = f [xi+k ... por tanto: f [xi .. por tanto ak = f [xi . .. . PN (x) = PN −1 (x) + aN (x − x0 )(x − x1 )... 3] = 2 e3 − 3e2 + 3e1 − 1 f [0.. xk ].. .. xi+k−1 ] f [xi .. x2 = 2... xi+k ] = xi+k − xi Teorema 12 Si denotamos por ak = f [x0 ... . 1... . ... ... 3] = e3 − e2 e2 − 2e + 1 f [0... 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
ak Πk−1 (x − xi ) i=0
donde los coeﬁcientes f [xi . 2. . xi+k ] = f [xi+k . .. por la unicidad del polinomio interpolador..... xi ] = bk De nuevo..... por tanto: ak−1 − ak
k−1 X j=0
Como P2 (2) debe ser igual a e2 . xk ]. xi+k ] − f [xi . P0 (x) = 1 P1 (x) = 1 + a1 x Como P1 (1) debe ser igual a e....(x − xN−1 ) A los coeﬁcientes ak los denotamos por ak = f [x0 . xk ] veriﬁcan f [xi+1 . cambiando el orden de los puntos.. entonces el polinomio de interpolación de Lagrange PN (x) viene dado por PN (x) =
Ejemplo 8 Sea f (x) = ex . si interpolamos f (x) en los puntos x0 = 0. xk de la siguiente forma: P0 (x) = a0 P1 (x) = P0 (x) + a1 (x − x0 ) P2 (x) = P1 (x) + a2 (x − x0 )(x − x1 ) .....
X.N) que a partir de los coeﬁcientes dados por el vector A(0 : N ) y el conjunto de puntos de interpolación.EQ.
PARAMETER(Nmax=1000) DIMENSION A(0:1000).’Puntos de Interpolacion repetidos’ STOP ENDIF PRINT *.B(0:Nmax) DO 1 K=0. aproximar f (x) por la parábola que pasa por los puntos (xn−1 .1) THEN PRINT *. 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 ). x1 = 3.F. π .X. x2 = 3 y x3 = 4 y.X(0:1000) PRINT *.N PRINT *.N) Parameter(Nmax=1000) DIMENSION A(0:*).25 0.-1 1 EVDIFNEWTON=EVDIFNEWTON*(X0X(K))+A(K) END
..A(K) PRINT *. 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 *.2 -0..15 -0.1 0.0.X0. Polinomio’ DO 3 K=0.N).
FUNCTION EVDIFNEWTON(A.X. x2 = 1 y x3 = 0.F. 3 y 4 utilizando las diferencias divididas de Newton. f (xn−1 )) .N DO 3 L=0. 3] :
0. Problema 30 (3 puntos) Dada una función f (x) y una secuencia de valores xn . Programa 10 Programa en Fortran 77 donde se deﬁnen las funciones IDIFNEWTON. x1 = −1. x2 = 1. f (xn−2 )) y (xn−3 .X(0:*).X.05 0 0 -0.1 -0.F(K). Ptos Interp.X(K). 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.EQ.F(0:1000). x3 = 2 utilizando las diferencias de Newton y evaluar el polinomio en x = 0 utilizando el algoritmo de Horner.5 x 3
3 Problema 27 (2 puntos) Interpolar la función f (x) = 10 x2 +1 en los puntos x0 = −2.F(0:*). 1.N B(K)=F(K) A(0)=F(0) DO 2 K=1.’ DO 1 K=0. . X(K) PRINT *.F(K) IF(IDIFNEWTON(A.N) DIMENSION A(0:*).N) END FUNCTION IDIFNEWTON(A. π y 3π utilizando las diferencias divididas 2 2 de Newton. devuelve el valor de la evaluación del polinomio de Lagrange en el punto X0.’Introducir Ptos Interpol. y la función EVDIFNEWTON(A. xK ] ). x1 .X(0:*) EVDIFNEWTON=A(N) DO 1 K=N-1. x1 = 1.’Coef. Expresar el polinomio tomando en primer lugar x0 = 0.’Test de Comprobacion’ DO 4 K=0.5 2 2.N READ *.’Introducir No.X.En la siguiente gráﬁca se muestra la diferencia ex − P3 (x) en el intervalo [0. devuelve el vector A(0 : N ) de coeﬁcientes de diferencias divididas que deﬁnen el polinomio de Lagrange (A(K) = f [x0 .05 -0.EVDIFNEWTON(A.X(K).X0. f (xn−3 )). en segundo lugar. (xn−2 .N-K IF (X(K+L). x0 = 4.’Introducir Valores de F()’ DO 2 K=0.N N=N-1 PRINT *.N PRINT *.’ READ *.5 1 1.X(L)) THEN IDIFNEWTON=1 RETURN ENDIF B(L)=(B(L+1)-B(L))/(X(K+L)-X(L)) CONTINUE A(K)=B(0) CONTINUE IDIFNEWTON=0 END
Problema 29 (3 puntos) Calcular el polinomio interpolador de Lagrange P3 (x) de la función f (x) = 2x en los puntos 0.
π ]. Asi. . 6
Comprobar que el polinomio esta bien construido. 2 y 3. 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).. lo que simpliﬁca el cálculo numérico.Implementación de funciones elementales Una vez deﬁnida una aritmética en precisión ﬁnita y las 4 operaciones básicas (suma. 8 × 10−9 4 (2 ∗ 5 + 1)! Por tanto.5. π ] (x) a las funciones trigonométri4 4 cas deﬁnidas sobre el intervalo [0. 1. por tanto. obtenemos ya la mejor aproximación posible de ex en el intervalo [0. donde m es un número entero y x0 ∈ [0. 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. los desarrollos de Taylor y los algoritmos de búsqueda de ceros (como √ vimos anteriormente para x). que calcula el polinomio interpolador a partir de los puntos y valores de interpolación. el error relativo es menor que 6. evaluar y mostrar P6 (x) y ex . Práctica 3 (Aproximación de ex .
i = 0.... 0. 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. Las técnicas elementales para deﬁnir estas funciones consisten en utilizar la interpolación polinómica.6 × 10−8 y.. Podemos deﬁnir en4 tonces las siguientes funciones: ½ cos[0.. π ] (x) si x ≤ 2 sen[0. . π ] (x) si x ≤ π 4 4 cos[0. π ] (x) = π 2 sen[0. denotemos por cos[0. de em .2π] (x) =
cos[0. π ] (x) si x ≤ 2 − cos[0. es decir que P6 (xi ) = exi para todos los xi . la función ln(x).π] (x) si x ≤ π −sen[0. multiplicación. π ] (x) si x ≤ 4 π π ( cos[0. 4 ] 2 − x) si x >
π 4 π 4 π 2 π 2 π 2 π 2
sen[0. es decir 7 puntos de interpolación. del mismo orden que la unidad de redondeo u en una aritmética de 32 bits. π ]. por un lado. como son: la raíz √ cuadrada x.).π] (2π − x) si x > π sen[0. π ] (π − x) si x > 2
cos[0. N 2 2N + 2 obtenemos que el error relativo veriﬁca que: µ ¶N+1 0 e 1 | ex − PN (x) | ≤ ex0 (N + 1)!2N 2 Para N = 6.
. si trabajamos con una aritmética de 32 bits. π ] (x) = 2 cos[0..π] (x) si x ≤ π cos[0. en el cálculo 0 de ex para x0 ∈ [0. y por otro. 2 horas) Crear una función en Fortran 77 que devuelva el valor de ex con x ∈ [0. utilizando algunas relaciones trigonométricas es suﬁciente deﬁnir las funciones cos(x) y sen(x) en el intervalo [0. que evalua el polinomio interpolador en un punto. la función ex . las funciones trigonométricas: sen(x) cos(x) y tan(x). Introducir por teclado un valor x. 1] en una aritmética de 32 bits. π ] ( 2 − x) si x > π 4 4 sen[0. Puesto que estas funciones son 2π periódicas.π] (x) = cos[0. − x2 x4 x2n + + . tomando un polinomio de grado N = 6.π] (2π − x) si x > π
El desarrollo en Serie de Taylor centrado en 0 del cos(x) es: cos(x) u Pn (x) = 1. 1]. + (−1)n 2 4! (2n)!
y el error máximo cometido por el desarrollo de Taylor en un punto x ∈ [0. y EVDIFNEWTON(. Utilizando como puntos de interpolación los asociados a los polinomios de Chebychev: µ µ ¶¶ 1 2i + 1 xi = 1 + cos π i = 0. π ] y a partir de ellas deﬁnir las funciones para 4 cualquier valor x (en radianes). π ] (π − x) si x > 2 sen[0. Por ejemplo. 1]. Efectivamente. la función xy .. a partir de estas operaciones. 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). cuya unidad de redondeo u es del orden de 10−8 .). Dado que 0 ex = em ex podemos descomponer el cálculo de ex en el cálculo.π] (x) = ½ sen[0. división).2π] (x) = 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. Utilizar x = 0. el valor 4 máximo del error se encuentra en x = π . Aproximación de la exponencial ex Un número real x siempre se puede expresar como x = m + x0 . etc. como tan(x) es creciente en [0.. Nota: Utilizar las funciones de an. resta. π ] (x) y sen[0. es necesario deﬁnir. El error relativo es | Pn (x) − cos(x) | (x)2n+1 ≤ tan(x) cos(x) (2n + 1)! Además. las funciones elementales que todos usamos.h IDIFNEWTON(.
utilizando únicamente su valor en el intervalo [0. A continuación. denominado pivote. Para obtener A0 y b0 se calcula. π ] 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. N 2 4 2N + 2 Dado que ln(1) = 0.. 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. que la matriz sea simétrica o deﬁnida positiva. pero que requieren. aunque no es de los más rápidos. = podemos reducir el cálculo del ln(x) al rango de valores 1 2 ≤ x ≤ 1. 1] 2 2 |x−1| ≤1 | ln(x) | Por tanto: | ln(x) − PN +1 (x) | 2N+1 ≤ | ln(x) | (N + 2) µ ¶N +1 1 1 4 2N
b0 N a0 N. .l
Problema 35 (2 puntos) Calcular el número de operaciones básicas (sumas.. algo que. Por tanto.. y u = (ui ) es el vector solución buscado.Problema 31 (3 puntos) Aproximar la función sen(x) en el intervalo [0. 1
l=k+1 a0 k.N k = N − 1. 1]. . . Una vez obtenidos la matriz A0 y el vector b0 . que para el intervalo [ 1 .. N. restas. 1] son: 2 µ µ ¶¶ 1 1 2i + 1 xi = + 1 + cos π i = 0.j ) es una matriz de N xN. en una aritmética de 32 bits tendríamos la mejor aproximación posible de la función ln(x). utilizando relaciones trigonométricas. Problema 34 (1 punto) ¿Cómo se puede obtener la función y x . tiene la gran ventaja de que se puede aplicar a todo tipo de matrices.. 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) ∼ . por tanto. para minimizar el error relativo añadiremos como punto interpolante xN +1 = 1. para que el sistema sea equivalente. b = (bi ) es un vector de tamaño N que determina los términos independientes. se intercambia la primera ﬁla de A con la ﬁla donde se encuentra el pivote. el valor máximo en valor absoluto de la primera columna de A. el número es mayor o n=1 2n igual que 1 y menor que 1. donde x e y son números reales. 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. como veremos en el futuro. a1 = 1 y para´ > 1 an = 0 ³P t an ó an = 1. 8
Problema 33 (3 puntos) Calcular los polinomios necesarios para interpolar las funciones trigonométricas cos(x) y sen(x) en el intervalo [0. π ].. A continuación. multiplicaciones y divisiones) necesarias para realizar un remonte como el presentado arriba en función de la dimensión N ..
Para N = 10 el error máximo es 3. Además se tiene que en el intervalo [ 1 .. 8
ANÁLISIS NUMÉRICO MATRICIAL I En esta primera sección dedicada a la resolución de sistemas de ecuaciones lineales. 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. estudiaremos los métodos directos clásicos para la resolución de un sistema de ecuaciones de la forma Au = b donde A = (ai. no ocurre con otros métodos más rápidos. es posible calcular las funciones sen(x) y cos(x) para cualquier x (en radianes). 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 . 973 6 × 10−8 . es decir. que tenga valores nulos de la diagonal hacia abajo. que es menor que la unidad de redondeo u y. por ejemplo.k
a0 ul k. Método de Gauss Este método. el cálculo de la solución u es inmediata. Se hace lo mismo para el vector
. 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.6931471806). y se hace lo mismo con el vector b. en primer lugar. Utilizaremos los puntos de interpolación generados por los polinomios de Chebychev. utilizando las funciones ex y ln(x)?
Problema 32 (2 puntos) Demostrar que. π ] en una aritmética de 32 bits.
N max) que resuelve un sistema. forma siguiente ⎛ ⎞ ⎛ 0 ⎝ 24 ⎠ → ⎝ 8
La descomposición de la matriz A lleva las siguientes fases: ⎛ ⎞ ⎛ ⎞ −2 −2 0 6 18 12 −−→ −− ⎝ 6 18 12 ⎠ pivoteo ⎝ −2 −2 0 ⎠ 3 11 7 3 11 7 ⎛ ⎞ ⎛ ⎞ 6 18 12 6 18 12 −− − − − − − − − − − − −→ ⎝ −2 −2 0 ⎠ ceros 1a colu mna ⎝ 0 4 4 ⎠ 3 11 7 0 2 1 ⎛ ⎞ ⎛ ⎞ 6 18 12 6 18 12 −− − − − − − − − − − − −→ ⎝ 0 4 4 ⎠ ceros 2a colu mna ⎝ 0 4 4 ⎠ 0 2 1 0 0 −1 el vector b se ha transformado de la ⎞ ⎛ ⎞ ⎛ ⎞ 24 24 24 0 ⎠→⎝ 8 ⎠→⎝ 8 ⎠ 8 −4 −8
Programa 11 Programa en fortran 77 que implementa el método de Gauss.N. u2 = −6.Nrow.EQ.’Una columna es toda cero’ ELSE PRINT *.N DO 2 L=1. En el caso en que se ha terminado correctamente.Nmax). u el vector solución. donde A es una matriz triangular inferior.B(Nmax).B(K) CONTINUE ELSE IF(M.NmaxGAUSS) THEN IGAUSS=3 PRINT *. un vector auxiliar N row.’Introducir matriz’ DO 2 K=1. 1
Problema 36 (2 puntos) Resolver por el método de Gauss el sistema¶ µ µ ¶ µ ¶ −1 2 x 3 = 2 −1 y 0 Problema 37 (3 puntos) Calcular el número de operaciones básicas necesarias para descomponer el sistema Au = b en el sistema A0 u = b0 utilizando el método de Gauss. y así sucesivamente hasta llegar a la mencionada matriz A0 .B(K) M=IGAUSS(A. y con ello habremos obtenido un sistema equivalente tal que la primera columna es cero de la diagonal hacia abajo. La función devuelve un valor entero M que indica si se ha terminado correctamente (M = 0) o incorrectamente (M = 1.Nrow(Nmax) PRINT *.N-1 XMax=ABS(A(Nrow(K).’DIMENSION DEL SISTEMA MAYOR DE LA PERMITIDA’ RETURN ENDIF DO 1 K=1.N PRINT *. la solución se devuelve en el propio vector b. Se considera el sistema ⎛ ⎞⎛ ⎞ ⎛ ⎞ −2 −2 0 u1 0 ⎝ 6 18 12 ⎠ ⎝ u2 ⎠ = ⎝ 24 ⎠ 3 11 7 8 u3
de la misma forma.Nrow(*) DIMENSION U(NmaxGAUSS) IF (N. teniendo en cuenta la siguiente relación:
M−1 X k=1
1 3 1 2 1 M − M + M 3 2 6
Problema 38 (2 puntos) Implementar en FORTRAN la funcion IDESCEN SO(A.’SOLUCION’ DO 1 K=1. ’Introducir Dimension’ READ *.
FUNCTION IGAUSS(A. 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. Volvemos ahora a hacer lo mismo para convertir la segunda columna cero de la diagonal para abajo.K)) M=K DO 3 L=K+1.B(*).Nrow.*). N.B. Se deﬁne una función IGAU SS que tiene como parámetros la matriz A.N.1)THEN PRINT *.N PRINT *.N READ *.’A(N.L) PRINT *. b. u3 = 8.N 1 Nrow(K)=K DO 2 K=1.B.A(K.b. la dimensión del sistema y la dimensión máxima admitida. Parameter(Nmax=1000) DIMENSION A(Nmax.
Problema 39 (2 puntos) Resolver por Gauss el siguiente sistema de ecuaciones: ⎛ ⎞⎛ ⎞ ⎛ 0 −1 2 u1 ⎝ −1 2 −1 ⎠ ⎝ u2 ⎠ = ⎝ 2 −1 0 u3
el método de ⎞ 1 0 ⎠ 1
Ejemplo 9 Ejemplo de descomposición según el método de Gauss.0) THEN PRINT *.EQ. u. La función devuelve 0 si termina correctamente y 1 en caso contrario.Nmax) PARAMETER (NmaxGAUSS=1000) DIMENSION A(Nmax.’Introducir vector’ DO 3 K=1.N READ *.GT .N
.N)=0’ ENDIF END
y el remonte da como solución u1 = 6.Nmax) IF(M. 2). el vector independiente b. b es el vector de términos independientes. Nota Importante: Las líneas de código tienen que ir todas numeradas y no pueden superar las 15 líneas.
.(2.**(-100))) THEN IGAUSS=2 RETURN ENDIF U(N)=B(Nrow(N))/A(Nrow(N).1 ⎜ b2.k k=1 Fin Para j Fin Para i
.**(-100))) THEN IGAUSS=1 RETURN ENDIF IF(K.. si la solución es perfecta entonces Au − b = 0.. Cuanto más pequeño sea ErrorSistema. Para estimar el error cometido al resolver el sistema utilizaremos la expresión siguiente.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. entonces A es simétrica y deﬁnida positiva.k bi.M)C*A(Nrow(K).IF(ABS(A(Nrow(L).. donde e es el vector e = Au − b : ErrorSistema = 1 X |ei | N |bi | + 1
Para j = i + 1..M) 5 CONTINUE B(Nrow(L))=B(Nrow(L))-C*B(Nrow(K)) 4 CONTINUE 2 CONTINUE IF(ABS(A(Nrow(N). En el denominador se añade 1 para evitar las posibles divisiones por 0. .i = bi. porque los errores de redondeo y de cálculo producen que esta estimación no sea exacta..
Método de Cholesky Este método sólo se puede aplicar a matrices simétricas y deﬁnidas positivas.i − i−1 b2 k=1 i.K)) M=L ENDIF 3 CONTINUE IF(XMax. ³ ´ P bi.1. el algoritmo para calcular B es el siguiente Para i =r .1 triangular inferior. Teorema 13 Sea A una matriz simétrica. mejor aproximada estará la solución del sistema. .M) THEN MP=Nrow(K) Nrow(K)=Nrow(M) Nrow(M)=MP ENDIF DO 4 L=K+1. bn. esto no suele suceder.2 0 0 b3.N C=C+A(Nrow(K).1 ⎜ B = ⎜ b3. bn.N) DO 6 K=N-1. un vector de términos independientes b y un vector solución u.K) 6 CONTINUE DO I=1.GT.1 ⎜ ⎝ . 0 .i = ai.N)). calculado utilizando alguna técnica numérica. 0 b2.3 . . (iii) Los determinantes de todos los menores principales de A son positivos. 0 .(2. 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.i − i−1 bj. El siguiente teorema da 3 posibles deﬁniciones equivalentes de una matriz deﬁnida positiva. N ³ ´ P 1 bj.LT.K)).2 . N 1.LT.. bn.n ⎞ ⎟ ⎟ ⎟ ⎟ ⎠
Problema 40 (2 puntos) Demostrar que si A = B · B t (B triangular inferior) y |B| 6= 0.K)/A(Nrow(K).L)*U(L) 7 CONTINUE U(K)=(B(Nrow(K))-C)/A(Nrow(K).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.N C=A(Nrow(L).2 b3.3 . 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. 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. bn.i aj.NE.K) DO 5 M=K.XMax) THEN XMax=ABS(A(Nrow(L).N A(Nrow(L). haremos lo siguiente: dada una matriz A. Ahora bien. .M)=A(Nrow(L).-1 C=0 DO 7 L=K+1. (ii) Todos los autovalores λ de A son positivos.
⎝ −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. Los 3 sistemas corresponden a matrices simétricas y deﬁnidas positivas.
Práctica 4 (Método de Cholesky. es muy sencillo resolver el sistema de ecuaciones Au = b. Estos archivos ejemplo sólo se pueden utilizar al compilar el programa en serdis bajo UNIX. a continuación. se almacena todo en una única matriz B.Nmax): Resuelve un sistema triangular superior.Nmax): Devuelve el error cometido al resolver el sistema dado por la expresión ErrorSistema de la sección anterior. Resolver los siguientes sistemas ejemplo: ⎞⎛ ⎞ ⎛ ⎞ ⎛ x 6 1 1 4 1.Nmax): Resuelve un sistema triangular inferior.Nmax): Calcula la descomposición de Cholesky de A y la devuelve en la matriz DB. .N.N. el programa no reconocerá el formato de los archivos. Si utilizamos Linux. Vamos a descomponer A en el producto de dos matrices triangulares de la forma siguiente: ⎛ ⎞ a1 b1 .dis.ulpgc. cN −1 aN ⎞ ⎞⎛ ⎛ 1 u1 .B. Efectivamente. ⎝ 1 5 6 ⎠ ⎝ y ⎠ = ⎝ 12 ⎠ 4 6 26 z 36 ⎛ ⎞⎛ ⎞ ⎛ ⎞ 1 1 4 x 6 2. Devuelve 0 si termina bien y −1 en caso contrario • FUNCTION ERROR_SISTEMA(A. ⎝ 1 5 6 ⎠ ⎝ y ⎠ = ⎝ 12 ⎠ 4 6 17 z 36 ⎛ ⎞⎛ ⎞ ⎛ ⎞ 2 −1 0 x 1 4.DZ. 0 0 l1 0 ⎜ m1 l2 0 ⎟ . Nota: El programa debe permitir introducir el sistema directamente por teclado o desde disco duro.h. Devuelve 0 si termina bien y −1 en caso contrario.DZ. tanto en doble como en simple precisión.B. Hay que resolver los sistemas ejemplo y calcular el ErrorSistema. donde DB es la matriz. donde DB es la matriz. utilizando las funciones deﬁnidas en an.DB. 0 1
Los vectores mi . para evitar tener que almacenar dos matrices.es. 0 ⎠ ⎝ 0 . ⎟ ⎟⎜ ⎜ ⎝ 0 . 100 y 500. simétrica. .El interés de descomponer una matriz A por el método de Cholesky es que. mN−1 lN 0 . . Problema 42 (2 puntos) Calcular el número de operaciones necesarias para resolver un sistema por el método de Cholesky. 0 ⎜ c1 a2 .N.
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).
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. 0 ⎟⎜ 0 1 . 6 horas) Implementar en Fortran 77 las siguientes funciones : • FUNCTION ICHOLESKY_FACTORIZACION (A. basta descomponer el sistema de la siguiente forma: Bz = b Bu = z
Ambos sistemas se resuelvan rápidamente haciendo un remonte y un descenso.DU. una para B y otra para B t . • FUNCTION IDESCENSO(DB. • FUNCTION ICHOLESKY(A. En este directorio hay tres ejemplos de sistemas de dimensión 10.N. • FUNCTION IREMONTE(DB. Devuelve 0 si termina bien y −1 en caso contrario. DZ es el término independiente y DU es el vector donde devuelve la solución. B es el término independiente y DZ es el vector donde devuelve la solución. N − 1
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.DU. bN −1 ⎠ 0 . . Nota: Normalmente.. y ui se calculan utilizando el esquema: l1 = a1 1 u1 = b1 l Para i = 2..DU.B.Nmax): Resuelve un sistema por el método de Cholesky y devuelve la solución en DU. Devuelve 0 si termina bien y −1 en caso contrario. . escribiendo en la parte triangular superior de B la parte correspondiente a B t . ⎝ 1 1 4 ⎠ ⎝ y ⎠ = ⎝ 12 ⎠ 4 4 17 z 36 ⎛ ⎞⎛ ⎞ ⎛ ⎞ 1 1 4 x 6 3. 0 ⎟ ⎜ ⎟= ⎝ 0 . uN −1 ⎠ 0 .N. li.
STATUS=’OLD’. Ndimension DO 1 K=1.Vector. STATUS=’OLD’.Vector) CHARACTER * (*) String DIMENSION Vector(*) OPEN(1.’=’. form=’UNFORMATTED’. form=’UNFORMATTED’.Ndimension.FILE=String.’=’. form=’UNFORMATTED’) 2 B=Ndimension WRITE(1) B DO 3 K=1.
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. ****** Función que lee un vector de disco de nombre String. ERR=1) GOTO 2 1 OPEN(1.STATUS=’OLD’.FILE=String.*) OPEN(1.Ndimension WRITE(1) Vector(K) 3 CONTINUE CLOSE(1) END
STATUS=’NEW’.Ndimension) CHARACTER * (*) String
***** Procedimiento que escribe en el ﬁchero String la matriz A ***** cuadrada de dimension Ndimension SUBROUTINE EscribirMatriz(String.’DIMENSION DE LA MATRIZ ’.Ndimension DO 4 L=1.A.Ndimension WRITE(1) A(K.FILE=String.A.FILE=String.L) 4 CONTINUE 3 CONTINUE CLOSE(1)
.ERR=1) GOTO 2 1 OPEN(1.Ndimension DO 2 L=1.
STATUS=’OLD’. form=’UNFORMATTED’) READ(1) B Ndimension=INT(B) PRINT *.Ndimension READ(1) A(K. FILE=String. form=’UNFORMATTED’) READ(1) A Ndimension=INT(A) PRINT *. STATUS=’NEW’.
DIMENSION Vector(*) OPEN(1.Ndimension DO 1 K=1. String. lo almacena ****** en la tabla A y devuelve la dimension de la matriz FUNCTION LeerMatriz(String.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.
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.Nmax) CHARACTER * (*) String DIMENSION A(Nmax.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 DEL VECTOR ’.FILE=String. lo almacena ****** en la tabla Vector y devuelve la dimension del vector FUNCTION LeerVector(String.Nmax) CHARACTER * (*) String DIMENSION A(Nmax. form=’UNFORMATTED’) 2 A=Ndimension WRITE(1) A 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.*) OPEN(1.String.
A. String deﬁne un string de caracteres de tamaño 10.Ndimension PRINT *. 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.V.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.h que se encuentra en el directorio de trabajo.1) de la función: ½ 2 x + y2 − 1 f (x. Consideremos la función f (x) = x3 . Se denomina orden de la aproximación a la potencia más pequeña que aparece en el término del error.String DIMENSION V(Nmax).K) CONTINUE PRINT * CONTINUE CALL EscribirMatriz(’matriz.K)=L+K PRINT *.
que está mucho más próximo al valor real.dat’. una fórmula de integración numérica es un procedimiento que permite aproximar el valor de la integral en un intervalo [a.Nmax) String=’vector. + (x − xi )N + .Nmax) Ndimension=LeerMatriz(’matriz.’N=’. Diferenciación Numérica La manera habitual de aproximar la derivada de una función f (x) en un punto xi consiste en utilizar el desarrollo de Taylor centrado en xi : f (x) = f (xi ) + f 0 (xi ) f N) (xi ) (x − xi ) + . diremos que el orden de aproximación es 1. 31 1. mientras que si xj < xi .A. obtenemos la siguiente expresión: f 0 (xi ) ≈ f (xj ) − f (xi ) + O (|xj − xi |) xj − xi
donde O (|xj − xi |) indica. La declaración CHARACT ER ∗ 10.1.N DO 2 k=1.Ndimension DO 6 L=1. en el cuerpo del programa.Ndimension A(L. la derivada se calcula hacia atrás.dat’ Ndimension=LeerVector(String.Ndimension DO 5 K=1. Por otro lado. básicamente.Ndimension DO 4 L=1.Ndimension) N=LeerVector(’vector3.h0 incluye el ﬁchero an.
Problema 47 (2 puntos) Calcular analítica y numéricamente la matriz gradiente en el punto (1.Ndimension PRINT *.A(L. b]. Si xj > xi .V) PRINT *. obtenemos la aproximación f 0 (1) ≈ 23 − 13 =7 2−1
Si tomamos ahora xj = 1. A(Nmax. obtenemos f 0 (1) ≈ 1.A(L.V(K) V(K)=2*V(K) CONTINUE CALL EscribirVector(’vector3. entonces la derivada se calcula hacia adelante. que el error cometido es una suma de potencias de |xj − xi | en la que la potencia más pequeña es 1. b] a partir de la evaluación de f (x) en algunos puntos incluidos en el intervalo [a.dat’.h’ PARAMETER(Nmax=100) CHARACTER * 10.Ndimension DO 1 K=1..1 − 1
Nota la declaración IN CLU DE 0an.V) PRINT *.dat’. truncamos el desarrollo de Taylor y despejamos.Ndimension.N PRINT *. 1! N!
Programa 13 Programa en Fortran 77 donde se describe un ejemplo de lectura/escritura de vectores y matrices. y) = x−y
. Ejemplo 10 Veremos en este ejemplo como. utilizando el valor de f (x) en otros puntos vecinos a xi . cuanto más próximo esté el punto xj al punto xi . 1) (utilizar h = 0. en este caso.V(K) CONTINUE Ndimension=3 DO 3 K=1.... INCLUDE ’an.13 − 13 = 3.K) CONTINUE PRINT * CONTINUE END
Si tomamos un punto xj 6= xi .Nmax) PRINT *. mejor será el valor aproximado de la derivada. Por lo tanto.dat’.
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 (xi − h). supondremos que la dimensión es 2.y) ∂y2
Si tomamos ahora xj = 0.y) Utilizaremos la siguiente nomenclatura: Fx = ∂F∂x .93 f 0 (1) ≈ = 3. Nótese que. la aproximación de las derivadas de una función de varias variables. lj). 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 demostrar que. y + l) = F + hFx + lFy + 1 (h2 Fxx + 2hlFxy + 2 ¡ ¢3 l2 Fyy ) + O( h2 + l2 2 )
xr − xl
f (xi )−f (xl ) xi −xl
aproxima la derivada segunda de f (x) en xi con un orden de aproximación de 1. calcular el polinomio de Lagrange que interpola a f (x) en esos 3 puntos. y − l) = F − hFx − lFy + 1 (h2 Fxx + 2hlFxy + 2 ¢3 ¡ l2 Fyy ) + O( h2 + l2 2 ) 7. 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. F (x − h. Para simpliﬁcar la exposición. y comprobar que da la misma fórmula que la presentada en el problema anterior. teniendo en cuenta que f (0) = 1. xi .y) ∂x∂y . se utilizan los desarrollos de Taylor siguientes en 2 variables. el orden de aproximación es 2. f (xi + h). Problema 50 (2 puntos) Calcular una aproximación de la derivada tercera f 000 (xi ) de una función f (x) en un punto xi . demostrar que la fórmula f (xi ) ≈ 2
00 f (xr )−f (xi ) xr −xi
∂ 2 F (x. Utilizaremos la notación Fi. Problema 51 (3 puntos) Dados 3 puntos. demostrar que la fórmula f (xi ) ≈
r )−f (x i )−f (x (xi − xl ) f (xxr −xi i ) + (xr − xi ) f (xxi −xl l )
Problema 52 (2 puntos) Considerar en el problema anterior que xl = xi − h. Problema 54 (2 puntos) Calcular una aproximación de la derivada primera y segunda de una función f (x) en x = 0. Deducir como queda la fórmula anterior para aproximar la derivada segunda. Problema 53 (3 puntos) Dados 3 puntos xl < xi < xr .
1. y xl = xi − h. la expresión anterior nos da f 0 (1) ≈ 23 − 03 =4 2
Diferenciación numérica en dimensiones superiores Estudiaremos. comprobamos que. f (xi − 2h). utilizando f (xi ).j ∼ = F (hi.y) . calcular el polinomio de Lagrange que interpola a f (x) en esos 3 puntos. en este apartado.Problema 48 (3 puntos) Dados 3 puntos distintos xl . y) = (Fx (x. calcular la derivada de ese polinomio en xi . F (x + h. f (4) = 9
aproxima la derivada de f 0 (xi ) con un orden de aproximación de 2. f (1) = 0. la precisión es mayor que con la fórmula anterior. F (x − h. F (x. Para discretizar las derivadas de una función F (x. y).y) . si tomamos xi = 1 y h = 1. y − l) = F − lFy +
h2 2 Fxx h2 2 Fxx
+ O(h3 ) + O(h3 )
l2 2 Fyy l2 2 Fyy
+ O(l3 ) + O(l3 )
6. y comprobar que da la misma fórmula que utilizando los desarrollos de Taylor. y). F (x + h. F (x + h. y) = Fxx (x.
. 01 0. y). Ejemplo 11 Veremos en este ejemplo como. más preciso es el valor de la derivada.13 − 0. y + l) = F + lFy + 4. que es el vector de derivadas parciales. xr . xr . y − l) = F + hFx − lFy + 1 (h2 Fxx − 2hlFxy + 2 ¡ ¢3 l2 Fyy ) + O( h2 + l2 2 ) 8. y el Laplaciano ∆F (x. utilizando la expresión anterior para aproximar la derivada de f (x) = x3 en x = 1. ∂x2
∂ 2 F (x. F (x − h. ∂y ∂ 2 F (x. si xr = xi + h. Fy (x. y)). Problema 49 (3 puntos) Dados 3 puntos distintos xl . y) = F − hFx + 3. (x. en este caso. Por ejemplo.
5. calcular la derivada segunda de ese polinomio en xi .2 que está más próximo al valor real que utilizando la primera fórmula.1 1. y) = F + hFx + 2. y xr = xi + h. y) + Fyy (x. cuanto mayor es el orden de una fórmula de aproximación. En general. Fy =
∂F (x. F (x. xi .
es decir: ∆F = Fi+1.j+1 + Fi+1.j 2h2 Fi+1. la función inicial en torno al punto (hi0 . un parámetro a elegir. 1 ) = 1 . Hablando en términos 3 de teoría de la señal. 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 .j−1 − 4Fi. hj0 ) = γ 2 + (1 − γ) 2 2h h Ahora bien. F (0. Por lo tanto.j−1 ) +γ 4h (Fi. F (0. hj0 ). − 1 ) = 2 2 2 2 2 2 2 1 1 1 1 1 1 1 2 .j−1 ) +γ 4h donde γ es. y) = (0. 0) = 1 .j−1 − 4Fi. F (− 1 . 2 ) = 0.j+1 − Fi. Problema 58 (2 puntos) Calcular una aproximación del gradiente de una función F (x. F ( 1 .Discretización del Laplaciano Para discretizar el operador ∆F en un entorno de 3 × 3 puntos. si queremos que ambos valores de ∆F coincidan. − 1 ) = 2 2 2 2 2 2 2 2 2 2 2 1 2. 0) = 0. 0) conociendo los siguientes valores: F (0. F ( 1 . 1 ) = − 1 . y) en el punto (x. F ( 1 . F (− 1 . 0) = 0.j−1 ) + 2h (Fi+1. F (− 2 .j + Fi. 0) = 1 .j−1 + Fi+1. La elección de dicho parámetro γ la haremos de forma que la discretización de ∆F respete lo máximo posible la invarianza por rotaciones de la función F (x.j+1 + Fi.j + Fi−1. 0) conociendo los siguientes valores: F (0. F (− 1 . hj0 ) = γ 2 + (1 − γ) 2 2h h Por lo tanto. debemos elegir γ = 2 .j+1 − Fi+1. cualquier promediado de las dos expresiones también es una discretización del laplaciano.j+1 − Fi−1. 1 ) = 1 .j+1 + Fi+1.j+1 + Fi−1.j ) + 2h (Fi+1. de donde deducimos. y).j h2
Problema 56 (2 puntos) Calcular una aproximación del laplaciano de una función F (x. − 1 ) = 1 . si rotamos 45 grados. lo será en particular para rotaciones de 45 grados.j+1 − Fi−1. − 2 ) = 0. consideremos una función tal que en un entorno de un punto (hi0 . y) = (0. 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. Para ello.j+1 + Fi−1. obtenemos como imagen: 1 1 0 1 0 0 0 0 0
Problema 57 (3 puntos) Demostrar que las máscaras √ ¢ √ ¢ ¡ ¡ − ¡√ 2 ¢ 0 2− 2− 2 ¢ ¡√ 1 Fx = −2¡ 2 √ 1 − ¢ 0 2¡ 2 √ 1 −¢ 4h − 2− 2 0 2− 2 √ ¢ √ ¢ ¡ ¡√ ¢ ¡ − 2− 2 −2 2 − 1 − 2− 2 1 0√ ¢ Fy = ¡ 0√ ¢ ¡√ 0 ¢ 4h ¡ 2− 2 2− 2 2 2−1
Si calculamos de nuevo ∆F en el mismo punto obtenemos: 1 2 ∆F (hi0 . F ( 2 . de nuevo. obtenemos la siguiente expresión para el gradiente: (Fi+1. F (0. Para simpliﬁcar. F (0. Problema 55 (3 puntos) Demostrar.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.j+1 + Fi.j =γ + 2h2 Fi+1.j + Fi−1.j )x = (1 − γ) (Fi. 1 ) = 1 .j − Fi−1.j + Fi.j−1 − 4Fi.j−1 − Fi−1. y) en el punto (x. que las siguientes expresiones son discretizaciones del laplaciano: ∆F ∆F = = Fi+1.j +(1 − γ) + h2 +O(h) donde γ es un parámetro libre a elegir.j+1 + Fi−1. por tanto. que γ = 2 − 2. utilizando el mismo √ argumento que para el ∆F . F (− 1 . 2 4 2 4 2 4 2 4 F ( 1 . − 1 ) = 1 . Discretización del gradiente Siguiendo el desarrollo de Taylor mostrado anteriormente. supondremos que l = h. F (− 2 .j−1 + Fi−1. el calculo de ∆F nos llevaría a convolucionar la imagen con la siguiente máscara: 1 h2
1 3 1 3 1 3 1 3 −8 3 1 3 1 3 1 3 1 3
dan lugar a una discretización del gradiente tal que su norma euclídea es invariante por rotaciones de 45 grados.j+1 + Fi−1. 2 ) = −1. 0) = − 1 . 2 2
. pueden utilizarse diferentes esquemas. 0) = 1 . − 1 ) = 1. F ( 1 . Teniendo en cuenta que la norma euclídea del gradiente es invariante por rotaciones.j−1 − 4Fi. utilizando el desarrollo de Taylor.
y Hn (x) = 2xHn−1 (x) − 2(n − 1)Hn−2 (x) para n ≥ 2. b].N los ceros del polinomio de x Legendre LN (x).8611363116 0. 1].5773502692 1 3 0. 0. Es decir. Z
cos(x)dx
Deﬁnición 3 Se denominan polinomios de Legendre Ln (x) a la familia de polinomios dada por L0 (x) = 1.. utilizando los ceros y pesos asociados a los polinomios de Legendre.5 y los pesos w0 = w1 = w2 = 2/3. b] utilizando la evaluación de f (x) en ciertos puntos de [a.5773502692 1. se utilizan los ceros de los denominados polinomios de Hermite... 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.3478548451
Cuando el intervalo [a. para cualquier polinomio P (x) de grado menor o igual que M.3478548451 0. utilizando las fórmulas de Legendre para n = 2 y n = 3: Z 1 ¡ 3 ¢ x − x4 dx
¿Cuál es el valor exacto de la integral?
f (x)dx ≈
wk f (xk )
donde xk representa los puntos de evaluación de f (x) y wk el peso de cada punto de evaluación. una fórmula de integración numérica se puede escribir como Z
Problema 59 (2 puntos) Aproximar el valor de la siguiente integral. la fórmula de integración numérica aproxima la integral de la siguiente forma: Z
f (x)e−x dx ≈
. es decir. . la fórmula es exacta. vamos a aproximar el valor de la integral de f (x) en [a. b].6251451549 − 0.8611363116 0.5. Demostración [Hu] Pg. cuál sería la fórmula de integración numérica de Legendre utilizando un sólo punto de interpolación. y para n = 2. Usar esta fórmula de integración para calcular númericamente la siguiente integral y compararla con el resultado análitico (exacto). para el intervalo [−1. 1].5555555556 0. −0.6251451549 −0. Si deﬁnimos wk = ˜ Z
Πi6=k (x − xi) ˜ dx Πi6=k (˜k − xi) x ˜
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.Integración Numérica Métodos de Cuadratura de Gauss Sea f (x) una función deﬁnida en un intervalo [a. En este caso. b] cualquiera. ∞). deﬁnidos como H0 (x) = 1. L1 (x) = x. En el caso [a. los puntos x0 = −0.. y dado un intervalo [a. ¿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. nLn (x) = (2n − 1)xLn−1 (x) − (n − 1)Ln−2 (x)
Problema 61 (2 puntos) Encontrar. a = −∞ o b = ∞.. b] es inﬁnito.3399810436 0.. b] = (−∞. 1]. H1 (x) = 2x.7745966692 0. x1 = 0 y x2 = 0.5555555556 4 0. Es decir Z b N X P (x)dx = wk P (xk )
a k=1
Problema 60 (2 puntos) Se consideran. 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. Deﬁnición 2 Una fórmula de integración numérica se denomina exacta de orden M si.3399810436 0. b]. 3. 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. hay que emplear otros métodos para aproximar las integrales. Estos puntos y estos pesos se utilizan para aproximar la integral de una función en [−1.8888888889 − 0.
se pueden utilizar los desarrollos a partir de los polinomios de Legendre. < xM+1 = b. lo cual resulta complejo para valores grandes de N. 211-213
¶ xk + xk+1 dx = 2 xk ¶ µ xk + xk+1 (xk+1 − xk ) = f 2 Z
. 707 106 781 0. ∞). ∞). L1 (x) = 1 − x. 1. En este caso. la fórmula de integración numérica aproxima: Z
donde a = x0 < x1 < . y Ln (x) = (2n − 1 − x)Ln−1 (x) − (n − 1)2 Ln−2 (x). para aumentar la precisión es necesario aumentar el grado de los polinomios. 585 786 438 0. deﬁnidos por L0 (x) = 1. Demostración [Hu] Pg..Teorema 15 Si xk son los ceros del polinomio de Hermite ˜ y deﬁnimos Z ∞ Πi6=k (x − xi) −x2 ˜ wk = ˜ e dx Πi6=k (xk − xi) ˜ −∞ entonces la fórmula de integración numérica generada por ˜ los puntos xk y los pesos wk es exacta hasta orden 2N − 1 ˜ para el intervalo (−∞. utilizando dos puntos de aproximación. el valor de la integral: Z ∞ 1 dx 1 + x2 −∞
M XZ
Para el intervalo (0.
Fórmulas de Integración Numérica Compuestas Con las fórmulas que hemos visto hasta ahora. 1. 886 226 925 5 0. Una alternativa consiste en dividir previamente la integral en subintegrales de la manera siguiente: Z
Problema 65 (2 puntos) Aproximar. 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 . 772 453 851 2 −0. 2 0. Demostración [Hu] Pg. ∞). 853 553 390 3 3. 414 213 562 0.
utilizando los polinomios de Hermite. 146 446 609 3
Problema 66 (2 puntos) Calcular de forma exacta la integral Z ∞ ¡ 3 ¢ x − x2 e−x dx
utilizando los polinomios de Laguerre.. A continuación se aproxima numéricamente cada una de las integrales Z xk+1 f (x)dx
Para ello. 707 106 781 0. 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. se utilizan los polinomios de Laguerre Ln (x). para n ≥ 2.
Problema 64 (2 puntos) Calcular de forma exacta la integral Z ∞ ¡ 3 ¢ 2 x − x2 e−x dx
Problema 67 (2 puntos) Calcular una fórmula de aproximación numérica de la integral siguiente: Z ∞ f (x)e−x dx
donde a es un número real cualquiera. 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.. 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.
la integral Z 1 ¡ 3 ¢ x − x4 dx
utilizando únicamente el valor de la función en los puntos: −1. 0. por el método de Simpson. Los parámetros de la función serán: Los límites del intervalo de integración en precisión real y el número de subintervalos en los que se dividirá el interFórmula de Simpson valo inicial. estudiaremos las técnicas de integración numérica sobre dominios Ω de dimensión superior a 1. 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.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 . 0 sin(x)dx = 2 Z xk+1 f (x)dx ≈ xk R1 x ¶ Z xk+1 µ 2.j
. Para realizar esta elección se utilizan técnicas de cuadratura. 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. tienen la ventaja de que pueden ser utilizadas cuando sólo conocemos la función a integrar en un conjunto equiespaciado de puntos. Es decir: 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.
donde debemos elegir los puntos (xi .j
Aunque estas fórmulas sean menos precisas que las deducidas a partir de los ceros de los polinomios de Legendre. Integración numérica en dimensiones superiores En esta sección. 1 y 1. 2 horas)
Ahora bien. Es decir. Nótese que. 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. yj )
Ω i. − 1 . teniendo en cuenta los resultados de la sección anterior sobre derivación numérica f 00 (xm ). supondremos que la dimensión es 2. pretendemos aproximar Z F (x. Es decir. se puede aproximar como f 00 (xm ) ≈ f (xk+1 ) − 2f (xm ) + f (xk ) ³ ´2
xk+1 −xk 2
Crear una función en fortran 77 donde se implemente el método de Simpson. 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. sarrollo en serie de Taylor centrado en el punto xm = xk +xk+1 . Para simpliﬁcar la exposición. y)dxdy
Por tanto. es decir. Es decir: 2 Rπ 1. sustituyendo este valor en la aproximación anterior obtenemos Z
f (x)dx ≈ f (xm )(xk+1 − xk )+ µ ¶ f (xk+1 ) − 2f (xm ) + f (xk ) xk+1 − xk + = 3 2 ³ ´ f (xk+1 ) + f (xk ) + 4f xk +xk+1 2 (xk+1 − xk ) = 6
Aproximaremos esta integral a través de la fórmula numérica: Z X F (x. 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. y)dxdy ≈ wij F (xi . La función a integrar se deﬁnirá aparte (como ³ ´ Z xk+1 en el caso del mérodo Müller). en este caso. yj ) y los pesos wij . la integración a partir de los ceros de los polinomios de Legendre no puede utilizarse. −∞ e−x dx = π = 1. 0 √1−x2 dx = 1 f 00 (xm ) 0 2 ≈ f (xm ) + f (xm )(x − xm ) + (x − xm ) dx = 2 xk µ ¶3 R∞ √ 2 f 00 (xm ) xk+1 − xk 3.
yk ) x e
1 3 y0 + y1 2 y0 + y2 y2 = e 2 y2 + y1 y3 = e 2 y1 = e
Problema 71 (2 puntos) Calcular de forma exacta la integral Z 1Z 1 x2 y 2 dxdy
x2 e tos. al igual que en dimensión 1. y [−1. 1] resultante de aplicar la integración numérica en una variable en los intervalos [−1. y2 ). y)e−x−y dxdy
x2 e x3 e x4 e
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. el área viene determinada por ¯⎞ ⎛¯ ¯ 1 1 1 ¯ ¯ ¯ 1 AREA(T ) = ABS ⎝¯ x0 x1 x2 ¯⎠ ¯ ¯ 2 ¯ y0 y1 y2 ¯
y. d] resultante de aplicar la integración numérica en una variable en los intervalos [a. que. los valores de los puntos y los pesos. y) = 3 3 T Integración sobre triángulos utilizando 3 puntos. Es decir P (x. el cálculo es un poco más complejo. podemos escribir: Z xm y n dxdy =
xm dx
y n dy
AREA(T ) el área del triángulo T. por los polinomios de Legendre. En función de los vértices. (x1 . la exactitud en dimensión 2 la podemos deducir a partir de la exactitud en dimensión 1. y) ∼ AREA(T ) =
donde w1 x1 e = w2 = w3 = = = = x0 + x1 2 x0 + x2 2 x2 + x1 2
wk F (ek . y) = ax + by + c. como hemos visto anteriormente. Deducir cual debe ser el punto (x0 . En este caso. b]x[c. Nótese que.
. d]. y)e−x
∞Z ∞ 0
F (x. µ ¶ Z ∼ F x0 + x1 + x2 . y)dxdy ≈ F (x0 . Consideremos un triángulo T de vértices (x0 . En el caso de que Ω sea un triángulo. x3 e
utilizando integración numérica. Problema 70 (2 puntos) Deducir la fórmula de integración numérica sobre un rectángulo [a. 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.
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. 1]. yk ) x e
F (x. d]. b]x[c. y) en 4 puntos. y0 )w0
sea exacta para polinomios de grado 1 en x e y. 1]x[−1. (x2 . Z F (x. y0 + y1 + y2 AREA(T ) F (x. 0). y0 ) y el peso w0 para que la fórmula de integración numérica: Z F (x. (1. en general. por tanto. 1]. 1). y) ∼ AREA(T ) =
donde w1 x1 e = w2 = w3 = = = = =
wk F (ek . 0) y (0. viene dada. y1 ). Denotaremos por
25 27 w4 = − 48 48 6x0 + 2x1 + 2x2 6y0 + 2y1 + 2y2 y1 = e 10 10 2x0 + 6x1 + 2x2 2y0 + 6y1 + 2y2 y2 = e 10 10 2x0 + 2x1 + 6x2 2y0 + 2y1 + 6y2 y3 = e 10 10 x0 + x1 + x2 y0 + y1 + y2 y4 = e 3 3
Problema 73 (2 puntos) Se considera el triángulo T de vértices (0. Un caso particularmente sencillo es cuando Ω es un rectángulo [a.donde m y n determinan el grado de los polinomios. y [c. en este caso. Problema 69 (3 puntos) Deducir la fórmula de integración numérica sobre el rectángulo [−1. De estas relaciones se puede deducir. y0 ). b]. también podemos extender los resultados al caso en que los intervalos sean inﬁnitos.
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. entonces. se podría deﬁnir su norma considerando la matriz como un vector de dimensión N xN . incluyendo técnicas iterativas de resolución de sistemas de ecuaciones y cálculo de autovalores.. (2. donde p es un número real positivo. 2). x2 ) que veriﬁcan que 1. en el espacio vectorial de los números reales.. Por ejemplo. Problema 79 (2 puntos) Demostrar que si A y B son dos matrices de dimensión N xN. k x k∞ < 1 Problema 78 (2 puntos) Tomar N = 2 y demostrar la siguiente desigualdad: k x k∞ ≤k x k2 ≤k x k1 Dada una matriz A de dimensión N xN . Entonces. Sin embargo. 0). 3 puntos y 4 puntos.
. 0) y (0. subordinada a la norma vectorial k . entonces.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 . k una norma vectorial. cuando trabajamos en varias dimensiones. se veriﬁca k AB k≤k A k · k B k
k x kp =
| xi |
Un caso particularmente interesante es p = 2. para cualquier norma de matrices subordinada a una norma vectorial. k x k1 < 1 2. la desigualdad es trivial. deﬁnida por k x k∞ = max | xi |
Problema 75 (4 puntos) Tomar N = 2 . p = 2. . y ∈ E. ANÁLISIS NUMÉRICO MATRICIAL II En esta sección veremos algunos aspectos más avanzados del análisis matricial. la norma ”natural” es el valor absoluto. x = (x1 . esta desigualdad es cierta por la propia deﬁnición de k A k. utilizando 1 punto. 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 . x2 . Básicamente. k x k2 < 1 3. esto es. que viene deﬁnida por ÃN X
1 !p
k Ax k kxk
La propiedad fundamental que veriﬁca una norma matricial deﬁnida de esta forma es la siguiente: Teorema 17 Sea A una matriz y k . k una norma vectorial.. y demostrar que la norma k x kp veriﬁca las propiedades de la deﬁnición de norma. lo que da lugar a la denominada norma inﬁnito. la desigualdad anterior es equivalente a k Ax k ≤k A k kxk Ahora bien. que corresponde a la norma euclídea. para cualquier vector x se veriﬁca que k Ax k≤k A k · k x k Demostración: Si x = 0. La deﬁnición más utilizada es la denominada norma p. Se deﬁne la norma de A. Si x 6= 0. una norma mide la magnitud o tamaño de un vector x. k como k A k= sup
donde Ω es el triángulo de vértices (0. puesto que k x k> 0. Otro caso interesante es aquél que se produce cuando hacemos tender p hacia inﬁnito. existen múltiples formas de deﬁnir una norma. xN ). Sin embargo. 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.
Normas de vectores y matrices Deﬁnición 4 Una norma k . Empezaremos recordando algunos conceptos relacionados con los autovalores..
tal que Ax = λx
xi . 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. por tanto. + (ηN λN ) k Ax k2 ≤ ρ(A) k x k2 para cualquier vector x. Puesto que A posee una base ortonormal de autovectores
Demostración: [La-Th] Pg.. k una norma vectorial. Deﬁnición 8 Se deﬁne el radio espectral de una matriz A como ρ(A) = max{| λi | : λi autovalor de A}
Y.. Teorema 19 Si los autovectores de una matriz A de dimensión N xN forman una base ortonormal de RN . al tomar el supremo en x. de la forma: x = η 1 x1 + η2 x2 + . obtenemos que
Deﬁnición 7 Se denomina polinomio característico P (λ) de la matriz A. + η N xN Al hacer Ax. denominado autovector. Teorema 21 Sea A una matriz cualquiera. En consecuencia.
Problema 82 (2 puntos) Calcular las normas 2. lo que demuestra que k A k2 ≤ ρ(A)
Problema 80 (1 punto) Demostrar que los autovalores de A son los ceros del polinomio característico P (λ). entonces todos sus autovalores son reales y. 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. entonces existe un autovector x tal que Ax = λx. al polinomio dado por el determinante P (λ) =| A − λI | que
Ax = η1 λ1 x1 + η2 λ2 x2 + . la desigualdad se mantiene. por tanto k Ax k k λx k = = |λ| ≤k A k kxk kxk Lo que demuestra el teorema. su producto escalar veriﬁca que (xi . Demostración: [La-Th] Pg.
Teorema 18 Sea A una matriz y k . Sea x un vector cualquiera. + (ηN )2 q 2 2 k Ax k2 = (η1 λ1 ) + . 1 e inﬁnito de la matriz µ ¶ 1 0 A= 1 1
.. 53. Vamos a demostrar la desigualdad k A k2 ≤ ρ(A) Dado que el teorema anterior determina la desigualdad en el otro sentido. en primer lugar. además. entonces k A k2 = ρ(A) Demostración: Recordamos. tendríamos la igualdad. se cumple q k x k2 = (η 1 )2 + . y puesto que los xi son autovectores. sus autovectores forman una base ortonormal de RN . Entonces k A k≥ ρ(A) Demostración: Si λ es un autovalor de A. + η N λN xN Como los autovectores son ortonormales. además. entonces p • k A k2 = ρ(t AA) P • k A k1 = maxj ( i | aij |) ´ ³P • k A k∞ = maxi | aij | j
(xi )k (xj )k =
0 si i 6= j 1 si i = j
donde (xi )k indica la coordenada k-ésima del vector xi . 73.. el vector x se podrá expresar como una combinación lineal de autovectores. xj ) =
Teorema 20 Si una matriz A de dimensión N xN es simétrica. y por tanto el resultado del Teorema.75.Deﬁnición 6 Un autovalor de A es un número λ real o complejo tal que existe un vector x.
entonces Limn→∞ k An k= 0 ⇐⇒ ρ(A) < 1 Demostración:[La-Th] Pg. Obviamente. 0.86. la perturbación de la solución del sistema puede llegar a ser del orden de 13. para la norma 2. a pesar de que la perturbación del sistema es del orden de 0. si los autovectores de una matriz A de dimensión N xN forman una base ortonormal de RN .1 ⎟ ⎜ 22.
. por tanto. 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. 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.01.1 ⎠ 30. al mismo tiempo.84.6.1). −12. de donde δu = A−1 δb y. 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. y 30.5. 3.9 ⎛
La solución de este sistema es (9. 1. Es decir. Problema 85 (2 puntos) Calcular el condicionamiento para la norma 2.6. 80. 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
Teorema 23 Si deﬁnimos χ(A) =k A k · k A−1 k entonces k δu k k δb k ≤ χ(A) kuk kbk
Demostración: Como A(u + δu) = b + δb y Au = b. Problema 84 (2 puntos) Demostrar que. 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.29.Problema 83 (2 puntos) Demostrar la siguiente igualdad: ρ(t AA) = ρ(A ·t A) Teorema 22 Sea A una matriz cualquiera.9 ⎟ ⎟=⎜ ⎟ ⎠ ⎝ 33. 4. ° ° kδuk ≤ °A−1 ° kδbk Por otro lado.1. 1. cuanto más pequeño sea χ(A). mejor comportamiento numérico tendrá la matriz A. también se cumple que kbk = kAuk ≤ kAk kuk de donde obtenemos que 1 kAk ≤ kuk kbk Así. Vamos a considerar ahora el mismo sistema. Como podemos observar. por tanto el condicionamiento sería χ(A) = 30.29 = 3029 0. 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. se obtiene que Aδu = δb.2. entonces.01
cuya solución es (1. multiplicando esta desigualdad con la anteriormente obtenida para kδuk . 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. Consideremos de forma genérica un sistema de ecuaciones de la forma Au = b y.
lo cual indica un condicionamiento bastante malo. 1). Para ilustrar de qué estamos hablando. concluimos la demostración del teorema.
707 107 ⎝ 0 3. π .0 0 −. sin α . la matriz R12 es 4 ⎞ ⎛ √ 1 1 − √2 0 2 1 1 √ 0 ⎠ R12 = ⎝ √
Para evitar tener que evaluar funciones trigonométricas. . 707 107 2. 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. que tienen la forma siguiente: ⎞ ⎛ 1 0 0 0 0 0 0 ⎜ 0 1 . 0 ⎟ Rpq (α) = ⎜ 0 . Concretamente. . 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.0 −. 1 . . cot(2α) = cos 2α (aqq − app ) = sin 2α 2apq
Ejemplo 16 Consideremos la matriz ⎛ ⎞ 2 −1 0 ⎝ −1 2 −1 ⎠ 0 −1 2
Para convertir en 0 el elemento a12 = −1 de la matriz.Método de Jacobi Este método se aplica a matrices reales y simétricas. 0 ⎟ ⎟ ⎜ . Además. sign(x) = 1 si x ≥ 0 y sign(x) = −1 4 4 si x < 0. 0 ⎟ ⎟ ⎜ ⎝ 0 . q apj sin α + aqj cos α j 6= p. dadas dos matrices A y R. 0 ⎟ ⎟ ⎜ ⎜ 0 . 707 107 ⎠ −. debemos elegir α tal que cot(2α) = (a22 − a11 ) =0 2a12
t y al hacer la operación R12 AR12 obtenemos ⎛ ⎞ 1. q
cot(2α) =
− tan(α) + sin(2α) 2 sin2 (α)
Utilizando las anteriores igualdades trigonométricas. 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 α ∈ − π . Este método intenta diagonalizar A realizando transformaciones del tipo R−1 AR. . Problema 86 (2 puntos) Sean las matrices A y R. se veriﬁca −1 que (Rpq (α)) =t Rpq (α). ⎟ ⎜ ⎜ 0 . q = apj sin α + aqj cos α j 6= p. Por tanto. basta con elegir α tal que pq (app − aqq ) sin 2α + apq cos 2α = 0 2
. Al realizar la operación A0 =t Rpq (α)ARpq (α).0
De donde α = − π . . Se basa en el hecho de que. y simpliﬁcar el algoritmo. haciendo 0 los elementos no diagonales mayores en módulo. sólo se ven afectadas las ﬁlas y columnas de índices p y q. . − sin α . cos α = 1 p 1 + tan2 (α) sin α = tan(α) cos α
donde los cosenos y senos están situados en las columnas y ﬁlas p y q. . si A es una matriz simétrica. q
El método de Jacobi se basa en ir modiﬁcando la matriz A mediante el procedimiento anterior. cos α . Al ser una matriz de rotación. la matriz A0 también es simétrica. 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. cos α . 1 0 ⎠ 0 0 0 0 0 0 1
Es decir. En el método de Jacobi se utilizan las denominadas matrices de rotación. que son costosas computacionalmente. se veriﬁca que los autovalores de A son los mismos que los autovalores de R−1 AR. Para anular un valor a0 . 707 107 −.
inicializamos B a la identidad antes de entrar en el bucle. j) FIN IF FIN PARA j A(p. . q) A(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. q) = A(q. j) = SI ∗ D + CO ∗ A(q.... Por tanto. j) A(j.Ahora bien. q) A(q.. + C ∗ C)
FIN IF CO = 1. p))/(2 ∗ A(p. DIM HACER IF ( j 6= p AND j 6= q) HACER D = A(p.. q) = A(q. . En primer lugar. j)) > R HACER R = ABS(A(i. q) = A(q. p) − T ∗ A(p. demostrar las igualdades siguientes: a0 pq a0 pp a0 qq a0 pj a0 qj = = = = = 0 app − tan(α)apq aqq + tan(α)apq apj cos α − aqj sin α apj sin α + aqj cos α
j 6= p. + T ∗ T ) SI = CO ∗ T PARA j = 1. Los parámetros de entrada son la matriz simétrica A.. · RM = D
donde D es una matriz diagonal que contiene los autovalores de A en la diagonal. haciendo cero el elemento no diagonal mayor en módulo de la matriz Ak−1 .. LOS AUTOVALORES SE ENCUENTRAN EN LA DIAGONAL DE A. Además los elementos no diagonales de A convergen hacia 0. + C ∗ C) ELSE T = −C + SQRT (1.. de tal forma que
−1 −1 RM · .. . RM . el número máximo de iteraciones N max y la tolerancia T OL para decidir cuándo son ceros los elementos no diagonales. la matriz B determina los autovectores. que se transforman de la siguiente manera: b0 ip b0 iq = cos(α)bip − sin(α)biq = sin(α)bip + cos(α)biq i = 1. DIM HACER PARA j = 1.Problema 89 (3 puntos) Dentro del método de Jacobi para el cálculo de autovalores. como Ri es una matriz de rotación del tipo Rpq (α). N i = 1.. j) A(j. 576-577.. para calcular la matriz B en el algoritmo anterior que calcula los autovalores.. Al utilizar el método de Jacobi. i − 1 HACER IF ABS(A(i.. vamos transformando la matriz A multiplicándola por una secuencia de matrices de rotación R1 . · R1 AR1 · . añadiremos en cada iteración las operaciones necesarios para ir obteniendo B. q
Veamos ahora la convergencia del método de Jacobi para el cálculo de autovalores. Numéricamente.. Despejando de la anterior igualdad obtenemos que AB = BD Si denotamos por bi el vector columna i de la matriz B. q)) IF C < 0 HACER T = −C − SQRT (1. A continuación... PARA n = 1.. entonces los elementos diagonales de la matriz Ak convergen (k → ∞) hacia los autovalores de la matriz A.. q j 6= p. q)) PARA i = 3. j) = CO ∗ D − SI ∗ A(q. y sea Ak la matriz transformada de Ak−1 . . j)) p=j q=i FIN IF FIN PARA j FIN PARA i IF R < T OL HACER PROCEDIMIENTO TERMINADO CORRECTAMENTE... Denotemos por B la matriz B = R1 · .. bi es el autovector de A asociado al autovalor dii .... · RM . N max HACER p=2 q=1 R = ABS(A(p. Algoritmo del Método de Jacobi para el Cálculo de autovalores. p) = A(p. q) + T ∗ A(p./SQRT (1. p) = 0 FIN PARA n PROCEDIMIENTO TERMINADO INCORRECTAMENTE NÚMERO DE ITERACIONES MÁXIMO EXCEDIDO Veamos ahora cómo podemos calcular los autovectores. q) − A(p. . su dimensión DIM. de la expresión anterior se obtiene que Abi = dii bi Es decir. . Teorema 24 Sea una matriz A simétrica. SALIR FIN IF C = (A(q. Sea A1 = A. Demostración: [La-The] Pg. en cada iteración haremos B = B · Ri . N
Problema 90 (3 puntos) Utilizar el método de Jacobi para aproximar los autovalores y autovectores de la siguiente matriz: ⎛ ⎞ 2 0 1 A=⎝ 0 1 0 ⎠ 1 0 1
.. p) = A(p.
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. ⎝ 2 ⎠ ↔ 2 − 2 ⎭ ⎩ 1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎛
10 y 100 del directorio de
• FUNCTION JACOBI(A.
⎞ 2. 626 28 . A = ⎝ 2 1 1 ⎠ −2 1 1 Resultado: ⎫ ⎧⎛ ⎧⎛ ⎞ ⎞⎫ 1 ⎬ ⎨ −2 ⎬ ⎨ ⎝ −1 ⎠ ↔ −2. 268 11
⎜ ⎜ ⎜ 16. 593 8 ⎜ −1.05 170 ⎟ ⎟ ⎜ 2. 6 ↔ ⎜ ⎜ ⎜ ⎝ ⎛
.0 −1. 998 09 −10. 400 89 ⎟ ⎠ ⎝ 1.B.N
max ERROR_V ECT ORES(A¯i . 659 86 ⎛ . 729 18 ⎟ ⎜ −. λ¯i ) x x
⎞ −. 741 24 ⎜ ⎜ −. Las matrices de dimensión la asignatura. entonces.Nota. Sea un vector u1 no ortogonal al subespacio engendrado por los autovectores del autovalor λmax .AUTOVALORES.N) : Devuelve la diferencia entre los vectores U y V. 176 3 1.N. 324 6 ⎟ ⎟ ⎜ 1. 165 3 ⎞ ⎛ −1.0 . ⎛ 0 1 6 0 0 0 ⎜ 1 0 2 7 0 0 ⎜ ⎜ 6 2 0 3 8 0 4. A = ⎜ ⎜ 0 7 3 0 4 9 ⎜ ⎝ 0 0 8 4 0 5 0 0 0 9 5 0 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠
⎧⎛ ⎧⎛ ⎞⎫ ⎞⎫ 1 ⎬ ⎨ −1 ⎬ ⎨ √ ⎝ 0 ⎠ ↔ 2. al comparar los resultados. x ¯ Devolver la expresión ERROR_AU T OV ECT ORES =
i=1. de dimensión N. 942 ↔ ⎜ ⎜ −1. Para no tener que buscar en cada paso el máximo de los elementos no-diagonales de Ak . 12 ↔ ⎜ −1. 853 75 ⎟ ⎟ ⎜ ⎜ . si deﬁnimos la secuencia un = A se veriﬁca que Limn→∞ sign
un−1 k un−1 k
¡¡ n n−1 ¢¢ u . utilizando la fórmula :
N 1 X ABS(U (i) − V (i)) ERROR_V ECT ORES = N i=1 ABS(U (i)) + 1. 215 57 ⎟ ⎠ ⎝ 1. ⎝ 1 ⎠ ↔ 2 ⎩ ⎭ 1 ⎛ ⎞ 2 −1 0 2. 901 1. 088 8 ⎟ ⎟ 5. 985 4 ⎟ ⎟ ⎜ ⎠ ⎝ 1. 521 32 ⎜ . 796 84 ⎟ ⎟ ⎜ ⎜ 3. 465 ↔ ⎜ ⎟ ⎜ −.0 −2.TOL. 938 54
Comprobar los resultados obtenidos con los siguientes ejemplos.Nmax): Para comprobar que los autovalores λi y su autovectores xi están bien estimados. La función devuelve 0 si termina bien y −1 en caso contrario. 826 8 ⎟ ⎟ ⎟ ⎜ ⎠ ⎠ ⎝ 1.u k un k= λmax
. los vectores A¯i y λi xi .V.Nmax. 11 ↔ ⎜ ⎟ ⎜ −.0001 y N iter = 1000: ⎛ ⎞ 2 2 −2 1. 896 18 ⎟ ⎟ ⎟ ⎜ ⎟ ⎟ ⎜ ⎟ − 12. • FUNCTION ERROR_AUTOVECTORES (A. 746 5 ⎜ .0 −1. que los autovectores están deﬁnidos módulo la multiplicación por una constante. utilizando la función ERROR_VECTORES().0 . tomando T OL = 0. 984 4 ⎜ −. el algoritmo de Jacobi se puede modiﬁcar haciendo cero el primer elemento apq que se encuentre que veriﬁque |apq | ≥ T OL. 855 2 ⎜ ⎝ 1. 575 36 ⎟ ⎟ ⎜ − 2.N.AUTOVECTORES. 562 8 ⎟ ⎟ ⎜ ⎜ 0. 241 2 ⎞ ⎞ ⎛ . ⎝ − 2 ⎠ ⎭ ⎩ ⎭ ⎩ 1 1 ⎧⎛ ⎞⎫ 1 ⎬ √ √ ⎨ √ 2.
Resultados: ⎛
3. B es una matriz donde se guardan los autovectores por columnas.Niter): Realiza el cálculo de los autovalores y autovectores de una matriz simétrica A por el método de Jacobi.0 . 06 ↔ ⎜ ⎜ −. A = ⎝ −1 2 −1 ⎠ 0 −1 2 Resultado:
Nota: Obsérvese. 363 27 . ⎝ −1 ⎠ ↔ ⎩ ⎭ ⎩ ⎭ 1 1 ⎞⎫ ⎧⎛ ⎨ 0 ⎬ 4. 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. comparar ¯ para cada autovalor λi .
+ µ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 + .¡¡ decir sign un . si aplicamos el método anterior a A−1 . + µN xN donde supondremos que x1 es un autovector asociado a λmax y que µ1 6= 0. |λmax | obtenemos que la secuencia tienden hacia 0. + µN |λλN | xN ° culo del autovalor de módulo menor λ ... un−1 = un = k un k
Para n = 1 la igualdad se cumple por la deﬁnición de u2 ... si λmax es negun = A−1 k un−1 k ativo. o (−1)n .. Por tanto. un−1 para n suﬁcientemente grande es 1n si λmax es positivo o (−1)n si λmax es negativo.. xM los autovectores asociados a λmax . un−1 es el signo del¢¢ producto ¡¡ escalar de un y un−1 . Problema 91 (3 puntos) Aplicar el método de la potencia para aproximar el autovalor máximo y el autovector asociado de las siguientes matrices. vamos a demostrar por inducción la siguiente igualdad: un+1 An u1 = kAn−1 u1 k
Limn→∞ kun k = ° ³ ° ´n−1 ³ ´n−1 ° ° λN °µ1 λmax x1 + ..... + µN |λmax | xN ° |λmax | ° ° °= = |λmax | ° ³ ´n−2 ³ ´n−2 ° ° °µ1 λmax x1 + . Demostración. y u1 no es ortogonal al espacio generado por los autovectores asociados a λmax . salvo si λi = λmax . como A posee una base de autovectores. En este caso. Sean x1 . entonces u1 se puede escribir como u1 = µ1 x1 + . un−1 k un k= λmin un =
.¡ ¡¡ ¢¢¢n Limn→∞ sign un . obtenemos que Limn→∞ sign ¡ ¡¡ ¢¢¢n Limn→∞ sign un .. un−1 < 0.. 1). 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 . Supongamos que se cumple para n−1... En primer lugar. + µN λn−1 xN max N °= = ° 1 n−2 °µ1 λmax x1 + . u ).. Por otro lado. ¡¡ ¢¢ Teorema 26 sign un . + µN λn−2 xN ° kAn−2 u1 k N ³ ´n−1 ³ ´n−1 Método de la potencia inversa λ x1 + .. dicho autovector tiene norma 1. µ ¶ 2 1 A = 0 1 µ ¶ −3 0 A = 1 1
An−1 u1 µ λn−1 x1 + . teniendo en |λmax | ° ° max min cuenta que Cuando hacemos tender n hacia inﬁnito. veriﬁca que Además ¡¡ ¢¢ 1 Limn→∞ sign un . un−1 = 1 si es ¡ n n−1 ¢ ¢¢ ¡ ¢ u . y demostrémoslo para n: un =A k un k
An−1 u1 kAn−2 u1 k kAn−1 u1 k kAn−2 u1 k
un+1 = A
An u1 kAn−1 u1 k
µ1 x1 + . realizando 3 iteraciones en el método.u k un k= λmax ¡ ¡¡ ¢¢¢n Por otro lado. + µN |λλN | xN ° |λmax | ° ° max = |λmax | y por tanto ¡¡ n n−1 ¢¢ u . + µM xM kµ1 x1 + . hasta calcular u4 y partiendo de u1 = (1. el término sign un . un−1 = −1 si un . que denotaremos por xi . si λmax es positivo... + µM xM k
Con lo que queda demostrado este primer resultado. Por la igualdad anteriormente demostrada obtenemos que
que es un autovector de λmax de norma 1.u ≥ 0 y sign un . un−1
un es un autovector de λmax k un k
Además... ... para n suﬁcientemente grande el signo de n n−1 λmax coincide con el signo del producto escalar (u .
Nótese que si el autovalor µ está calculado con mucha precisión. 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.. 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.h IGAUSS(). 1. Por ejemplo. 1) y tomando como norma kuk = maxi |ui |.. 1. se evita calcular directamente A−1 . podemos aplicar el método de la potencia inversa anterior. Para evitar esto. haremos 1 = 10−11 A = A − µId J = IGAU SS(A0 . utilizando el método de Jacobi. 0). podemos perturbar ligeramente el valor de µ para que IGAUSS() no dé problemas. Algorítmicamente. 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. y se obtiene un resolviendo el sistema Aun = un−1 k un−1 k
¡¡ n n−1 ¢¢ u . y. Para autovalores que se encuentren entre λmin y λmax . al hacer iteraciones de la forma un = M un−1 + c se obtenga que un converge hacia u.0) T HEN = ∗ 10. 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. 1). y como el determinante de una matriz es el producto de sus autovalores. entonces el autovalor más pequeño de A0 está muy próximo a 0. Si hacemos iteraciones del esquema anterior a partir de la aproximación inicial u1 = (0.u
Para ello. realizando 2 iteraciones del método a partir de u1 = (1.N E.
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.. la solución exacta del sistema es u = (1. y. quedaría como sigue: Si µ es el autovalor que estamos tratando. 1). calcular dos iteraciones del método de la potencia inversa partiendo de u1 = (1. 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 . por tanto. y1 . entonces se obtiene que el autovalor menor de A0 es justamente λ.Limn→∞ sign
un es un autovector de λmin k un k En los casos prácticos. obtenemos que x2 y2 z2 = = = 1+0 1 = 2 2 0+0 =0 2 1+0 1 = 2 2
Problema 94 (3 puntos) Calcular el autovalor y autovector más cercano a 2 de la matriz ⎛ ⎞ 0 −1 0 ⎝ 0 3 −1 ⎠ 0 0 −1
. µ = µ(1 + ) GOT O 1 EN DIF
Buscar la solución de este sistema es equivalente a buscar un vector u = (x. donde µ es uno de los elementos diagonales de la matriz que resulta de aplicar el método de Jacobi. 1).. 0.. la solución del sistema original. . 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.) IF (J. si consideramos la matriz A0 = A − µI.
Este método consiste en tomar MGS cGS = (D + L) = (D + L)
(−U ) b
A efectos prácticos. el vector candidato inicial uf. L es la matriz triangular inferior que corresponde a la parte de A situada por debajo de la diagonal... obtenemos que ⎛ ⎞ ⎛ ⎞ ⎛ ⎞ 0.Tolf ) que devuelve el autovalor máximo de una matriz y su autovector por el método de la potencia... − aNN −1 un−1 + bN 1 2 N −1 = aN N =
Si hacemos un barrido para el cálculo de la solución de arriba hacia abajo.. y en el caso de Jacobi se actualiza sólo al ﬁnal. −aN1 un − aN 2 un . Por tanto.Nf.84 .5 . Los parámetros son la matriz Af.Nﬁter.. − a2N un−1 + b2 1 3 N = a22 ..... Tomar como norma kuk = i ABS(ui ) (28 líneas de código como máximo). − a2N un−1 + b2 1 3 N = a22 . básicamente... entonces u veriﬁca que u = Mu + c
De la misma forma. 1).. 0. − a1N un−1 + b1 2 N a11 −a21 un−1 − a23 un−1 . u17 = ⎝ . y U es la matriz triangular superior que corresponde a la parte de A situada por encima de la diagonal. y vamos actualizando las componentes del vector aproximación según las vamos calculando.**120 si P termina no correctamente. En este caso. y Tolf la tolerancia. Ejemplo 18 Vamos a aplicar el método de Gauss-Seidel al sistema del ejemplo anterior.Nf )
Escribir siguientes: que
= 1 = 0 = 1
.. la aplicación de este método no requiere el cálculo directo de la matriz inversa (D + L)−1 .. donde D es la matriz diagonal que corresponde a la parte diagonal de A. después de haber calculado todas las componentes por separado.. puesto que el paso de una iteración a otra puede hacerse de la siguiente forma: un 1 −a12 un−1 − . Nf la dimensión real. Nfmax. 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. las sucesivas iteraciones se van aproximando a la solución u = (1. − aNN −1 un −1 + bN 1 2 N = aN N =
Método de Jacobi un 2 Este método consiste en tomar MJ cJ = D−1 (−L − U ) = D−1 b un N
Es el que se ha utilizado en el ejemplo anterior.Como puede observarse. 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. Esta función devuelve el valor 2. u8 = ⎝ .uf. −aN1 un−1 − aN 2 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).5 . 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. − a1N un−1 + b1 2 N a11 −a21 un − a23 un−1 . 1.73 ⎠ . y la función AUTOVALOR_MAXIMO(Af..25 ⎠ .96 ⎠ 0. Todas se basan en descomponer A de la forma A = L + D + U.84 . 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 − . obtenemos el método de Gauss-Seidel.. 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.Nfmax. Nﬁter número máximo de iteraciones.98 u3 = ⎝ 0. la dimensión para coger memoria.vf. es decir 2x − y −x + 2y − z −y + 2z
Problema 95 (3 puntos) en Fortran las funciones SIGNO_PRODUCTO_ESCALAR(uf.
en el caso de matrices tridiagonales. 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.. en general. matrices con todos los elementos nulos salvo la diagonal principal y sus codiagonales. 999 ⎛
= w (D + wL)
((1 − w)D − wU ) b
Estas nuevas matrices permiten realizar un promediado entre el resultado obtenido por Gauss-Seidel y el
. 0). en este caso. 976 56 ⎠ .. 999 u2 = ⎝ . 686 . u3 = ⎝ . Se toman. 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. y −1 c = (D + U ) b. Demostración [La-Th]. un problema difícil. 988 28
Problema 97 (2 puntos) ⎛ Calcular 1 0 de autovectores de la matriz ⎝ 0 2 1 0
una base ortogonal ⎞ 1 0 ⎠ 1
Como puede observarse de la expresión anterior. 342 ⎠ . 921 . obtenemos que ⎞ ⎛ ⎞ ⎛ ⎞ . 999 ⎠ . qué diferencias de implementación habría con respecto al caso anterior.... Sin embargo. 0.5 . −aN1 un . Ejemplo 19 Vamos aplicar el método de relajación al sistema del ejemplo anterior.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. 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. 0) y tomando w = wopt = 1. 976 56 u3 = ⎝ . u8 = ⎝ . ρ(MJ ) = √2 y wopt = 1. de la forma siguiente: un 1 un 2 un N −a12 un−1 − . Indicar..358-362. 17. 0)
Problema 99 (1 punto) Una variante del método de Gauss-Seidel consiste en tomar M = (D + U )−1 (−L). − a2N uN + b2 = w + (1 − w)un−1 2 a22 . u8 = ⎝ . Pg.. − a1N un−1 + b1 2 N + (1 − w)un−1 1 a11 n−1 n −a21 u1 . 0. es decir. − 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.17.. 625 . 802 ⎠ . obtenemos que ⎛ ⎞ ⎛ ⎞ . Mw cw = (D + wL)
1 + yn−1 + (1 − w)xn−1 2 xn + zn−1 = w + (1 − w)yn−1 2 1 + yn = w + (1 − w)zn−1 2 = w
Si hacemos iteraciones del esquema anterior a partir de la aproximación inicial u1 = (0.. 785 .. 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.
1 En este caso.. el valor de wopt se encuentra siempre entre 1 y 2. Teorema 28 Si A es una matriz tridiagonal y ρ(MJ ) < 1... 25 ⎠ . 585 . el siguiente resultado muestra la forma de calcular el valor óptimo de w.. 0. en este caso.
Demostración: El resultado es inmediato a partir del hecho de que una matriz M n converge hacia 0 cuando n → ∞ si y sólo si ρ(M ) < 1 Teorema 31 Si en el método de relajación w ∈ (0. 2). un .1 − aN−1.346-347.
entonces el método de Jacobi asociado al sistema Au = b converge para cualquier aproximación inicial. 0 − aNN−1.2 a ⎜ − aN −1.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. 0. 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. [La-The] Pg.N − aN. su determinante es el producto de los elementos diagonales. 2).1 aN. como el determinante de una matriz es el producto de sus autovalores.N . entonces los métodos iterativos convergen.
con la desigualdad estricta en al menos una ﬁla o columna. por tanto. Calcular previamente el parámetro de relajación óptimo.
| ajj |≥
| aij | ∀j. esto es ⎛ ⎞ 2 −1 0 ⎝ −1 2 −1 ⎠ 0 −1 2 satisface las hipótesis del Teorema anterior. Ejemplo 20 La matriz del sistema ejemplo tratado anteriormente. . − N. . entonces |1 − w| ≥ 1 y. Por tanto. obtenemos que.N −1 aN. Mw posee al menos un autovalor de módulo mayor o igual que uno. que kMJ k < 1 para la norma 1 o inﬁnito. observamos que la matriz MJ puede expresarse como: ⎛ 0 − a12 − a13 . Demostración. 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. − a1N a11 a11 a11 a21 a23 ⎜ − a22 0 − a22 . Demostración: En primer lugar. | aii |>
. Además. se tiene.N ⎝ aN−1. Entonces en = M n−1 e1 Demostración: La solución del sistema satisface que u = M u + c.
|Mw | =
|(1 − w)D − wU | (1 − w)N Πi aii = |(D + wL)| Πi aii
Por lo tanto.N−1 . Teorema 33 Si A es una matriz irreducible y se veriﬁca que X | aii |≥ | aij | ∀i. una matriz es irreducible si el cambio de cualquier valor del vector b del sistema Au = b afecta a todos los elementos del vector u. si w ∈ / (0. ⎜ aN −1. − a2N ⎜ a22 ⎜ . por las condiciones del teorema. 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. Demostración: En primer lugar. Teorema 29 Se considera el esquema iterativo un = M un−1 + c. Restando esta igualdad de la igualdad un = M un−1 + c.N −1 0 aN. el teorema se concluye teniendo en cuenta que cualquier norma de una matriz es siempre mayor o igual que su radio espectral. en consecuencia. teniendo en cuenta que el determinante del producto de dos matrices es el producto de sus determinantes y que el determinante de la matriz inversa es el inverso del determinante. Teorema 32 Si una matriz A veriﬁca que X | aij | ∀i. . .2 a − aN. / entonces ρ(Mw ) ≥ 1. observamos que las matrices D +Lw y (1−w)D −wU son matrices triangulares y.N−1 −1.partiendo de u1 = (0. obtenemos que
Problema 100 (3 puntos) Calcular 3 iteraciones del método de relajación para resolver el sistema ⎛ ⎞⎛ ⎞ ⎛ ⎞ 1 −1 0 x −1 ⎝ −1 2 0 ⎠ ⎝ y ⎠ = ⎝ 3 ⎠ 0 −1 3 z 1
o | ajj |>
| aij | ∀j.
. Problema 102 (2 puntos) Demostrar que.Nf.Problema 101 (2 puntos) Escribir en Fortran la función siguiente: CONDICIONAMIENTO(Af. Los parámetros son la matriz Af. uN ) = 0 donde f (u) = (f1 (u).. uN ) = 0 . reales o complejas.. Si el método no converge devuelve −1..*120 si termina mal porque Jacobi da un error o se produce una división por cero. si u1 y c son combinaciones lineales de autovectores de M correspondientes a autovalores de módulo menor que 1. donde z = x + yi. el parámetro de relajación w. Se supondrá implementada la función JACOBI(A. se escribe como N ecuaciones del tipo f1 (u1 . . 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..Niter) que devuelve 0 si termina bien y 1 si termina mal. Por ejemplo. Nﬁter número máximo de iteraciones. que por defecto se tomará 0. un sistema no lineal de ecuaciones de dimensión N. Nfmax. entonces el determinante de A es cero. de un polinomio de grado 2 dado por P2 (z) = az 2 + bz + c... ⎝ −1 2 −1 ⎠ ⎝ y ⎠ = ⎝ 0 ⎠ 0 −1 2 z 1 ⎞⎛ ⎞ ⎛ ⎞ 1 3 3 x 7 3. el número máximo de iteraciones N max. y por tanto el sistema asociado a A no tiene solución... la dimensión para coger memoria. ⎝ −1 2 0 ⎠ ⎝ y ⎠ = ⎝ 3 ⎠ 0 −1 3 z 1 ⎛ ⎛ ⎞⎛ ⎞ ⎛ ⎞
donde u0 es una aproximación de la solución de f (u) = 0.. Los parámetros de la función serán: la matriz A. un vector u donde se almacenará la solución.. calcular las raíces. si una matriz A veriﬁca que por ﬁlas o columnas su suma es siempre igual a 0... muchas veces nos encontramos con sistemas no lineales de ecuaciones. entonces el método converge.. 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. .. 2 horas) Desarrollar una función en Fortran 77 donde se implemente el método de relajación. uN ) = 0 f2 (u1 .Nﬁter) que devuelve el condicionamiento de una matriz utilizando el método de Jacobi para calcular los autovalores.
Método de Newton-Raphson para sistemas de ecuaciones no lineales En las aplicaciones reales. y la tolerancia T OL para evaluar la diferencia entre un y un−1 . el vector b. es equivalente a resolver el sistema ax2 + bx − ay 2 + c = 0 2ayx + by = 0 que es un sistema no lineal de ecuaciones.. aunque el radio espectral de M sea mayor que 1. fN (u)) es una función de <N → <N . 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. En general.. y u = (u1 . 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. Los sistemas ejemplos del directorio de la asignatura.. 1).TOL. y Tolf la tolerancia (21 líneas de instrucciones como máximo)..
4. y que inicialmente será el vector aproximación inicial. Probar el método para los sistemas 1 −1 0 x −1 1. La función CONDICIONAMIENTO devuelve 2. Nf la dimensión real.Nfmax. . 1.Nmax. fN (u1 . La función devolverá el número de iteraciones necesarias para alcanzar la solución.TOLf. . uN ).. .N. . Estos ejemplos tienen siempre como solución el vector (1. Comparar la diferencia en la velocidad de convergencia entre el método de Gauss-Seidel y el Método de relajación. ⎝ 3 1 3 ⎠ ⎝ y ⎠ = ⎝ 7 ⎠ 3 3 1 z 7 ⎛
Ejemplo 21 Consideremos el siguiente sistema no lineal de ecuaciones: x2 − y 2 + 1 = 0 2xy = 0 La matriz gradiente de esta función viene dada por µ ¶ 2x −2y ∇f (un ) = 2y 2x
Práctica 7 (Método de relajación... f2 (u).
Interpolación de Hermite En ocasiones. 1). Por tanto. como indica el ejemplo siguiente. tenemos que resolver el sistema µ 1 ¶µ ¶ µ ¡ ¢2 ¡ ¢2 ¶ 1 z1 −3 − 3 +1 2 2 4 4 =− 3 1 6 z2 2 2 16 ¡ ¢ 9 cuya solución es − 13 . El 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. utilizando el método de NewtonRaphson.. Interpolación por splines cúbicos Uno de los problemas básicos del polinomio interpolador de Lagrange.N ..Si partimos de u1 = (1. f 0 (a) f N ) (a) (x − a) + . 1).. − 40 . y) = (1. 1).
El error de interpolación viene dado por la fórmula f (x) − PN (x) = f N +1) (ξ) (x − a)N +1 (N + 1)!
donde ξ es un valor intermedio entre x y a.
Problema 108 (3 puntos) Calcular los polinomios base de Hermite que corresponden a tomar como puntos de interpolación x0 = −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. 1). En el caso general. 1). las raíces complejas o reales de un polinomio de grado 3. para valores grandes de N. + (x − a)N 1! N!
PN (x) = f (a) +
que ya es una buena aproximación de la solución exacta dada por el vector (0. sino también el valor de sus derivadas.j (x). u2 viene dado por µ ¶ µ 3 ¶ µ 1 ¶ 1 −4 2 4 u = + = 3 1 −1 4 4 Para calcular u3 . y los resultados obtenidos por la interpolación pueden no ser muy satisfactorios. y el orden de derivación M = 1. el tema de interpolación de funciones se dividió en dos partes. se utilizan los denominados polinomios base de Hermite Hi. 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. z) = (1. es que. resulta de interés interpolar no sólo el valor de la función en ciertos puntos {xi }i=0. 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. calcular u2 y u3 utilizando el método de Newton-Raphson para aproximar un cero del sistema no lineal. 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. el polinomio interpolador de Hermite se deﬁne como: P (x) =
N M X X ∂j f i=0 j=0
Problema 106 (2 puntos) Se considera el sistema no lineal (x − 1)y = 0 (y − 2)x = 0
(xi )Hi. siendo ésta la segunda parte. y. 1. Un ejemplo clásico de ello es el desarrollo de Taylor de una función en un punto a. x1 = 1.. 1). que son polinomios de grado menor o igual que (N + 1)(M + 1) − 1 dados por las siguientes condiciones: ½ ∂ l Hi. 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. para obtener u2 tenemos que resolver µ ¶µ ¶ µ ¶ 2 −2 z1 −1 = 2 2 −2 z2 ¡ 3 1¢ que tiene por solución − 4 . − 4 .j (x)
A partir de u1 = (1.. Por motivos de coordinación entre los programas teórico y práctico de la asignatura.N .. u3 viene dado por 40 µ 1 ¶ µ 13 ¶ µ ¶ 3 − 40 − 40 3 4 u = + = 3 9 39
Tomemos como aproximación inicial u1 = (1. los polinomios de grado N pueden tener un carácter fuertemente oscilante.
. Problema 105 (2 puntos) Plantear el algoritmo necesario para calcular. aproximamos f (x) por un polinomio de grado N .
cN−1. 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... Si hay N +1 puntos. entonces ai di bi = f (xi ) i = 0.Ejemplo 22 El polinomio base de Lagrange centrado en 0 sobre los puntos xi = −5.. xi+1 ]..5
ci+1 − ci 3hi
i De la Condición P3 (xi+1 ) = f (xi+1 ). −1. Por razones técnicas. c0 . de la condición obtiene que
i−1 ∂P3 ∂x (xi ). 0.. como veremos posteriormente. despejando. d0 .25 0 -5 -2.
bi = 3di−1 h2 + 2ci−1 hi−1 + bi−1 i−1 y. Para añadir estas dos ecuaciones hay dos procedimientos estándares.. deﬁniendo un polinomio distinto para cada intervalo [xi ... aN −1 . el número de polinomios es N.... N − 1 3hi ai+1 − ai hi (2ci + ci+1 ) = − hi 3
i = 0. N − 1.. se obtiene de forma inmediata que ai = f (xi ).. cuando se trabaja con muchos puntos de interpolación. 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. . 4. Nótese que... .5 x 5
1. De la condición
i+1 ∂ 2 P3 ∂x2 (xi+1 )
3. . 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. La técnica más conocida es la interpolación por splines cúbicos.... es decir: a0 . ..
2. ... Para deﬁnir estos polinomios. obtenemos que
0 2.. −4. . despejando todo en función de ci . que son polinomios de grado 3. . N − 2 ∂x i+1 ∂ 2 P3 (xi+1 ) i = 0. ..
i Teorema 34 Si una familia de polinomios P3 (x) = di (x− 3 2 xi ) + ci (x − xi ) + bi (x − xi ) + ai . 5].75
i ∂ 2 P3 ∂x2 (xi+1 ).. El primero consiste simplemente en ﬁjar c0 = cN = 0. .. satisface
Nótese que esta última relación determina un sistema de ecuaciones donde las incógnitas son las variables ci . N. Ejemplo 23 (x2 − 1)(x2 − 4)(x2 − 9)(x2 − 16)(x2 − 25) = −14400
las condiciones anteriores. se imponen las siguientes condiciones:
i P3 (xi ) = f (xi ) i = 0. dN−1. b0 . −2.. . N − 1
Despejando. vamos a utilizar también los valores aN y cN . Para completar dicho sistema.. N − 1 i P3 (xi+1 ) = f (xi+1 ) i = 0. 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 ∂P3 (xi+1 ) = ∂x i ∂ 2 P3 (xi+1 ) = 2 ∂x
i+1 ∂P3 (xi+1 ) i = 0... 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 . 3. .. cN ) y N −1 ecuaciones.. −3. 1.
i Demostración De la condición P3 (xi ) = f (xi ). bN−1 .. para deﬁnir los polinomios.. .. lo que signiﬁca que c0 cN
0 ∂ 2 P3 (x0 ) = 0 2 ∂x N ∂ 2 P3 −1 (xN ) = 0 ∂x2
2ci+1 = 6di hi + 2ci de donde. xi+1 ].... N ci+1 − ci = i = 0. N − 2 ∂x2
= Vamos a introducir la notación hi xi+1 − xi . Dicho sistema tiene N +1 incognitas (c0 . tenemos que buscar 4N valores. . i = 0. 2. hay que añadir una ecuación que involucre a c0 y otra ecuación que involucre a cN . obtenemos que bi = ai+1 − ai − di h2 − ci hi = i hi ai+1 − ai hi (2ci + ci+1 ) − hi 3
i ∂P3 ∂x (xi )
Finalmente. Por tanto. se suele interpolar la función utilizando polinomios a trozos..
6 + 0 = 2− = 0. cuales son los puntos de unión entre los tres polinomios.25 1 0.75 2.933 3
y 5 4 3 2 1 0 0 -1 x -2 -3 -4 0.8 c2 Los valores bi y di vienen dados por d0 d1 d2 = = = −2. Los bj y dj se calculan directamente a partir de las relaciones mostradas en el teorema anterior.75 3
. se resuelve un sistema de ecuaciones tridiagonal para el cálculo de los ci .2 = 1.933 (x − 2) + 2.733 3 −4.133 (x − 2)
Por lo tanto.El segundo procedimiento se utiliza cuando utilizamos los valores de f 0 (a) y f 0 (b).667 (x − 1)3 − 2.2 = 1.25 2 1.75 1.467 (x − 1) + 1
de donde salen las ecuaciones f 0 (a) = f 0 (b) = a1 − a0 h0 − (2c0 + c1 ) h0 3 aN − aN−1 hN −1 − (2cN −1 + cN ) hN −1 3
P0 (x) = −0.25 0.5
2. parece.5
0. tal y como se muestra en la gráﬁca de la derivada segunda siguiente:
cuya solución es ¶ µ ¶ µ c1 −2.75 2 2. sobre la gráﬁca de la derivada segunda los puntos de unión se detectan en los lugares donde encontramos un pico.4 + 2. imponemos que
0 ∂P3
N−1 ∂P3
(x0 ) = f 0 (a)
−2. 3].8 + 2. sabiendo que f (0) = 0. tampoco sobre la derivada se aprecian los puntos de unión de los polinomios. P1 (x). En este caso. Debemos deﬁnir 3 polinomios distintos que corresponden a los intervalos [0. no es posible distinguir geométricamente. f (3) = 2.5 0.25 0 -0.667 3 0. [1. 2].75 1 1. f (1) = 1.8 = −1 − = −0.25 1. y P2 (x). el trazado de una única función.5 0. f (2) = 0.733x
P2 (x) = −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. los polinomios son P1 (x) = 1. Sin embargo. para calcular los splines cúbicos es necesario. por las condiciones sobre las derivadas que hemos impuesto. y 3.25 2.75 0.5 2.2 − 0 = −0.8 (x − 2) + 0.25 1.75 -1 -1. Veamos ahora gráﬁcamente el perﬁl de la derivada de los polinomios P0 (x). tomar ai = f (xi ).8 = −0. A continuación.75 3 x
El sistema que debemos resolver para calcular los ci es µ 4 1 1 4 ¶µ c1 c2 ¶ µ −6 9 ¶
como puede observarse. siguiendo con el resultado del teorema anterior.733 3 2.5 1. en primer lugar.5 1. al trazar la curva. a simple vista.25 0.5 2. y [2.
2. 2. En este caso hi = 1.133 3 = 1−
(xN ) = f 0 (b)
Por tanto.733x3 + 1.25
0. 3] : Como puede observarse.2 (x − 1)2 −0.25 -0. 1. Ejemplo 24 Vamos a calcular los polinomios interpoladores utilizando splines cúbicos al interpolar la función f (x) en los puntos x = 0.75 2
2. tomando c0 = c3 = 0.5 -0.2 = . 2.5
1.25 2. − 2.75 1
1.467 3 5. Es decir. 1].
π . la función f (x) = sin(x) en los puntos x = −π. Problema 110 (2 puntos) Calcular la función que interpola.75 2 2. 2.5 0. y la interpolación utilizando la función sin c(x) (línea a trozos). 1. la interpolación por splines cúbicos es la menos oscilante.2 0 0 0.5
0.75 3 x
Ejemplo 26 Consideremos la función f (x). 2 2 La interpolación trigonométricos a través de polinomios
cuya gráﬁca es
La base de la transformada de Fourier discreta es la utilización de los polinomios trigonométricos dados por la expresión P k (x) = eikx
0. cuando el número de puntos de interpolación aumenta.6
Esta función tiene la propiedad de que en x sin c(0) = 1. − π .5
0 0 0. 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). interviene
1.75 1 1.. 2 La interpolación a través de la función seno cardinal Una base de funciones interpolantes muy utilizada en la teoría de Fourier es la base formada a partir de la función seno cardinal. deﬁnida por sin c(x) = cuya gráﬁca es sin(x) x
1. y 3.4 0. 1 y 2.2 1 0. utilizando la función sin c(x)..8 1.25 0.25 0 -50 -25 0 25 x 50
2 1. π. deﬁnida en los puntos x = 0.5 1 1. f (2) = 0. los polinomios de la interpolación por splines cúbicos (línea sólida).75
Ejemplo 27 Vamos a comparar gráﬁcamente el resultado de interpolar la función del ejemplo anterior utilizando la interpolación de Lagrange normal.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. f (3) = 2. tal que f (0) = 0. de 0. N dada por la función ¢ ¡ sin(π x − i ) e ¡x a ¢ f (xi ) f (x) = π a −i i=M
= 0. 0. y para cualquier entero i distinto sin c(πi) = 0.5 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.8 0.4 1. Dada una función f (x).25 2. Por otro lado.6 0. f (1) = 1.5 2 2. .5 2. la diferencia entre los diferentes tipos de interpolación también lo hace.25 1. su función polante en los puntos xi = a · i para i = M. la interpolación por splines cúbicos y la interpolación a través de la función seno cardinal.
a través de una función.. y por tanto ! Z π Ã N X ∂E ilx (c−N ..N . pretendemos aproximar f (x) como f (x) ≈
P3 (x) =
1 2 2 + cos(x) − cos(3x) 2 π 3π
ck eikx
La siguiente gráﬁca muestra la aproximación entre f (x) y P3 (x):
donde ck son coeﬁcientes. cN ). observamos que.Dada una función f (x). el polinomio trigonométrico interpolador es
Demostración En primer lugar. π]. . sin exigir que la función aproximante pase exactamente por ese conjunto de puntos. cN ) = cl e f (x) − eikx dx = 0 ∂ck −π
l=−N
Problema 111 (3 puntos) Calcular el polinomio trigonométrico.. 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. π]. que interpola la función f (x) = |x| en el intervalo [−π. observamos que. π ] 2 2 f (x) = 0 si x ∈ [− π . dada la forma cuadrática del funcional E(c−N . cN ) con respecto a cualquier ck son cero. deﬁnida en el intervalo [−π...25 0 -2. . ... 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 .. Dado un conjunto de valores {(xi . en general complejos. π ] / 2 2
(axi + b − yi )2
Vamos a calcular el polinomio trigonométrico interpolante para N = 3. Aproximación por mínimos cuadrados La aproximación mínimo cuadrática aproxima. teniendo en cuenta que ½ Z π 2π si l = −k eilx eikx dx = 0 si l 6= k −π
Ejemplo 28 Consideremos la función ½ 1 si x ∈ [− π . Por otro lado. yi )}i=1.. tal que la función de error cuadrático E(a.. tomando N = 2.75
0.. las derivadas parciales de E(c−N ... dada la forma cuadrática que tiene el funcional. 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. la aproximación mínimo cuadrática lineal consiste en buscar la recta y = ax + b. b) = sea mínima.5 x
vienen dados por ck = Rπ f (x)e−ikx dx −π 2π
Demostración En primer lugar.5
0. debe poseer un
. cN ) =
−π k=−N
0. .25 0 1... en un mínimo.5 -1.25 2. un conjunto de valores de forma global. éste debe poseer mínimos.
Masson . ”Programación en fortran 77” Anaya. Trae una buena selección de problemas. "Análisis Numérico". da una visión general sobre los últimos avances en Análisis Numérico. las derivadas parciales son cero. 1988. haciendo especial énfasis en la precisión de los algoritmos numéricos y en la propagación de errores. [Ci] Ciarlet P. 1990. 1 Méthodes directes y Vol. Su mayor virtud es el rigor matemático con el que se tratan los temas y una cuidada presentación. muestra las últimas tendencias en cuanto a la enseñanza de los conceptos básicos del Análisis Numérico. Inc.W. Masson. y por tanto X ∂E (axi + b − yi ) xi = 0 (a. 1989. iterativos y métodos tipo gradiente. [La-Th] Lascaux P. Los algoritmos están muy bien descritos a través de un seudocódigo. 1994. 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. Uno de los libros clásicos más conocidos en Análisis Numérico. 1966. la buena presentación de los temas elegidos la hacen de interés. sin pretender ser exhaustiva. con múltiples ejemplos y una descripción de los algoritmos bien diseñada. Vol. Destaca por el rigor matemático en su exposición. así como el cálculo de autovalores y vectores propios. incluyendo la resolución de sistemas a través de métodos directos. 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.G. Esta obra.E. b) = 2 ∂b i=1
Esto da lugar a un sistema lineal de ecuaciones cuyas incógnitas son a y b. Además. Esta obra. muy reciente.. Esta obra.mínimo. ”Análisis Numérico”. "Accuracy and Stability of Numerical Algorithms". "Analyse numérique matricielle appliquée à l’art de l’ingénieur. [St] Stewart G. [Hi] Higham N. Cheney W. Faires D. ”Afternotes on Numerical Analysis” SIAM. ”Analysis of Numerical Methods”. The Benjamin/Cummings Publishing Company. Addison-Wesley Iberoamericana.
.. trata en profundidad todos los tópicos relacionados con el Análisis Numérico Matricial. presenta una cuidada selección de temas básicos en Análisis Numérico. Con un exquisito rigor se abordan los temas básicos del Análisis Numérico Matricial y métodos de optimización. Keller H.. Contiene todos los tópicos habituales con una descripción muy completa y detallada. destaca por una exposición simple y al mismo tiempo clara. b) = 2 ∂a i=1
X ∂E (axi + b − yi ) = 0 (a. ”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. y cuya resolución lleva al resultado establecido en el teorema. [Ki-Ch] Kincaid D.
BIBLIOGRAFÍA BÁSICA [Bu-Fa] Burden R. F. 1993. 1996. Excelente libro de base para un curso de Métodos Numéricos. b).E. Grupo Editorial Iberoamérica 1985. John Wiley and Sons. 2 Méthodes itératives ". Théodor R. Esta obra es un clásico del Cálculo Numérico. en un mínimo del funcional E(a.E. [Hu] Hultquist P. 1996 Esta obra.. ”Introduction à l’analyse numérique matricielle et à l’optimisation”. [Bo] Borse G. dividida en dos volúmenes. sin pretender ser tan exhaustiva como otras obras de carácter más general. SIAM..
APÉNDICE C: Algunos fallos comunes en Fortran 1. Este comando es de utilidad para salvaguardar la información de directorios y ﬁcheros de miradas ajenas. No poner EN D al ﬁnal del programa principal o de la función.name :q! sale del vi sin guardar cambios.f chmod cambia los permisos de lectura. A continuación. Finalmente.f practica1. :w f ichero. escritura y ejecución de un ﬁchero. 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. y el modo edición.
.f ls (dir) visualiza contenido de un directorio >ls /users/p701/fortran77 cp (copy) copia un ﬁchero en otro. chown cambia el propietario de un ﬁchero. aparece el comando UNIX. donde se ejecutan comandos. Puede ejecutarse en dos modos. >cp /users/p701/fortran77/programas/prog1. Hacer > man chown para mirar las opciones.f . 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. cd (cd) cambia el directorio activo >cd /users/p701/fortran77 more (type) visualiza el contenido de un ﬁchero >more /users/p701/fortran77/programas/prog1.name escribe el ﬁchero actual en el ﬁchero f ichero. Hacer > man chmod para mirar las opciones.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.
Intercambio entre modo comando y modo edición ESC pasa de modo edición a modo comando i pasa de modo comando a modo edición A pasa a modo edición y pone el cursor al ﬁnal de la línea O inserta una nueva línea. pasa a modo edición y pone el cursor al principio de la nueva línea Manejo de Ficheros (en modo comando) :w escribe en disco el ﬁchero :wq escribe en disco el ﬁchero y sale del vi :e f ichero.name edita el ﬁchero f ichero. un comentario y un ejemplo. que es donde se escribe normalmente el texto. >mv prog1. El modo comando (el que está por defecto al entrar en vi). entre paréntesis. su equivalente en MS-DOS (si existe). rm (del) borra un ﬁchero >del prog1.f man (help) suministra ayuda sobre un comando > man ls logout se termina la sesión y se sale del sistema >logout ps visualiza los números de procesos que están abiertos que corresponden al usuario alumno >ps -u alumno kill interrumpe la ejecución de un proceso de número N proceso >kill -9 N proceso mkdir (mkdir) crea un directorio >mkdir practica1 rmdir (rmdir) borra un directorio >rmdir practica1 mv (move) cambia de nombre o ubicación un archivo.APÉNDICE A: Resumen de los comandos de UNIX En este breve resumen seguiremos el siguiente esquema. En primer lugar.
Solución: Escribir 1. 3. Utilizar variables enteras como ﬂotantes o al revés. No respetar los tipos en los pasos de parámetros de las funciones. 10.
. los programas pueden fallar por errores de redondeo en los cálculos. si en lugar de utilizar la directiva -rn utilizamos -dn aumentaremos la precisión de las variables declaradas DOUBLE PRECISION.f -o prueba” donde n es 8 ó 16. Sugerencia: Aunque no sea necesario. Escribir números como 1/2 ó 10 ∗ ∗20 en precisión entera. si hacemos ”f77 -rn prueba. 11. Solución: Poner una declaración de PARAMETER al principio de la función y con ella asignar las memorias de forma estática. Siempre hay que buscar que el número de anidamientos sea mínimo. Análogamente. 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. Las sentencias GOT O pueden diﬁcultar el seguimiento del ﬂujo del programa y sólo hay que utilizarlas cuando sean indispensables. Utilizar vectores sin declararlos con la sentencia DIMENSION. Exceso de sentencias GOT O. declarar los tipos de todas las variables que se utilicen al principio del programa o función. 8. 6. Por ejemplo. 9. y si utilizamos -in las variables enteras./2. No pasar la dimensión de un vector como parámetro de una función. Fortran da la posibilidad de cambiar el número de bits utilizados para almacenar las variables en el momento de la compilacion. 7. aumentaremos la precisión de la aritmética para las variables reales. A veces. No poner ningún comentario en los programas.2. 5.ó 10. 4. ∗ ∗20. Anidar excesivamente los programas.
Más de este usuarioPHS 100-1A21Daniel Nuñez HernandezPHS 100-1A21ARTICULADOS- INTRODUCCION 2014Daniel Nuñez HernandezARTICULADOS- INTRODUCCION 2014Angular Contact Ball Bearings, Single RowDaniel Nuñez HernandezAngular Contact Ball Bearings, Single RowDialnet IngenieriaInversaDeUnReductorDeTornilloSinfincoron 4271876 (1)Daniel Nuñez HernandezDialnet IngenieriaInversaDeUnReductorDeTornilloSinfincoron 4271876 (1)Maxon Motor 118392Daniel Nuñez HernandezMaxon Motor 118392ReductorDaniel Nuñez HernandezReductorMSDS-HPC Daniel Nuñez HernandezMSDS-HPC Isometrica TuberiasDaniel Nuñez HernandezIsometrica TuberiasTD11-03 Grafico de espacios maximos entre guias - (Macoga).pdfDaniel Nuñez HernandezTD11-03 Grafico de espacios maximos entre guias - (Macoga).pdfClave cengel tranfe.docxDaniel Nuñez HernandezClave cengel tranfe.docx68917550-Guia-de-Valvulas.pdfDaniel Nuñez Hernandez68917550-Guia-de-Valvulas.pdfFlujo compresible.pdfDaniel Nuñez HernandezFlujo compresible.pdf26_evaluacion.pdfDaniel Nuñez Hernandez26_evaluacion.pdfAqui estaba la papa.pdfDaniel Nuñez HernandezAqui estaba la papa.pdfManualCompostajeparaMunicipios.pdfDaniel Nuñez HernandezManualCompostajeparaMunicipios.pdfTest A la caldera Una vez finalizado el calorimetro.pdfDaniel Nuñez HernandezTest A la caldera Una vez finalizado el calorimetro.pdfManualCompostajeparaMunicipios.pdfDaniel Nuñez HernandezManualCompostajeparaMunicipios.pdfManual-de-Compostaje-Casero-MMA.pdfDaniel Nuñez HernandezManual-de-Compostaje-Casero-MMA.pdfỘģỘģ庨ħd - Analisis analitico.pdfDaniel Nuñez HernandezỘģỘģ庨ħd - Analisis analitico.pdfEnsamblaje1Daniel Nuñez HernandezEnsamblaje1Whats New EnterpriseDaniel Nuñez HernandezWhats New EnterpriseSW2012_SP4.0_fixed_SPRsDaniel Nuñez HernandezSW2012_SP4.0_fixed_SPRsNCh0003-1961Daniel Nuñez HernandezNCh0003-1961NCh0003-1961Daniel Nuñez HernandezNCh0003-1961ClavesDaniel Nuñez HernandezClaves

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