Universidad de Almería

Departamento de Estadística y Matemática Aplicada

Cálculo Numérico para I.T.I.G.

Apuntes

José Antonio Rodríguez Lallena

Almería, 24 de febrero de 2011

Índice general

Introducción 5

1. ANÁLISIS DE ERRORES 71.1. Números máquina . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71.2. Redondeo, desbordamiento y cancelación . . . . . . . . . . . . . . . . . . . . 10

1.2.1. Errores de redondeo . . . . . . . . . . . . . . . . . . . . . . . . . . . 111.2.2. Errores de cancelación . . . . . . . . . . . . . . . . . . . . . . . . . . 141.2.3. Errores de desbordamiento . . . . . . . . . . . . . . . . . . . . . . . . 15

1.3. Error absoluto y error relativo . . . . . . . . . . . . . . . . . . . . . . . . . . 161.4. Ejercicios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

2. RESOLUCIÓN NUMÉRICADE ECUACIONES Y SISTEMAS NO LINEALES 232.1. Introducción a los números complejos . . . . . . . . . . . . . . . . . . . . . . 232.2. Soluciones de ecuaciones y su multiplicidad . . . . . . . . . . . . . . . . . . . 26

2.2.1. Ecuaciones polinómicas, raíces y su multiplicidad . . . . . . . . . . . 262.2.2. Raíces y su multiplicidad en ecuaciones no polinómicas . . . . . . . . 29

2.3. Métodos de bisección, de la regula falsi yde la secante . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 312.3.1. Método de bisección . . . . . . . . . . . . . . . . . . . . . . . . . . . 312.3.2. Método de la regula falsi . . . . . . . . . . . . . . . . . . . . . . . . . 342.3.3. Método de la secante . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

2.4. Método de Newton-Raphson . . . . . . . . . . . . . . . . . . . . . . . . . . . 372.5. Método de Newton-Raphson para sistemas

de ecuaciones no lineales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 392.6. Ejercicios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

3. RESOLUCIÓN NUMÉRICADE SISTEMAS DE ECUACIONES LINEALES 493.1. Preliminares de Álgebra Lineal . . . . . . . . . . . . . . . . . . . . . . . . . 493.2. Normas matriciales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 543.3. Condicionamiento de una matriz . . . . . . . . . . . . . . . . . . . . . . . . . 563.4. Métodos iterativos de resolución de sistemas lineales . . . . . . . . . . . . . . 583.5. Métodos de Jacobi y de Gauss-Seidel . . . . . . . . . . . . . . . . . . . . . . 59

3.5.1. Método de Jacobi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

3

3.5.2. Método de Gauss-Seidel . . . . . . . . . . . . . . . . . . . . . . . . . 613.6. Método de las potencias . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 633.7. Ejercicios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

4. INTERPOLACIÓN DE FUNCIONES 734.1. Problemas de interpolación y de aproximación de funciones . . . . . . . . . . 734.2. Interpolación polinomial: la fórmula de Newton . . . . . . . . . . . . . . . . 754.3. Interpolación mediante funciones spline . . . . . . . . . . . . . . . . . . . . . 784.4. Ejercicios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82

5. INTEGRACIÓN NUMÉRICA 875.1. Fórmulas simples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 875.2. Fórmulas compuestas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88

5.2.1. Método de los trapecios . . . . . . . . . . . . . . . . . . . . . . . . . 895.2.2. Método de Simpson . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89

5.3. Estudio del error de las fórmulas compuestas . . . . . . . . . . . . . . . . . . 905.3.1. Error de la fórmula de los trapecios . . . . . . . . . . . . . . . . . . . 905.3.2. Error de la fórmula de Simpson . . . . . . . . . . . . . . . . . . . . . 92

5.4. Ejercicios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94

Bibliografía 99

4

Introducción

El objetivo principal de estos apuntes es facilitar el trabajo de los alumnos (que, porejemplo, no necesitarán tomar nota en clase de las de�niciones, resultados, etc.) y, en parti-cular, el estudio de la asignatura.

Los apuntes desarrollan de manera concisa y concreta la teoría de la asignatura, que seilustra con algunos ejemplos; contienen un número más que su�ciente de ejercicios propuestospara que se resuelvan �a mano, con calculadora o con ordenador� en clase o en el estudiopersonal de cada alumno. Las relaciones de ejercicios terminan con algunas cuestiones teóricasque se han preguntado en exámenes de la asignatura.

Los apuntes no contienen la parte de la asignatura que se desarrolla con el programaMathematica, pero sí tratan de algunas cuestiones relacionadas, y proponen ejercicios quedeben ser resueltos con la ayuda de dicho programa.1 Una introducción sencilla a Mathemat-

ica puede encontrarse en Ramírez González el al. [10]. El manual de Mathematica se incluyeen la ayuda del programa, pero también puede ser útil consultar los manuales de Wolfram[18, 19]. Unas guías más completas acerca de Mathematica son las de Trott [14, 15, 16, 17].

¾Y qué es el Cálculo Numérico? El Cálculo Numérico �denominado también Análisis

Numérico� es una rama relativamente joven de la Matemática2 que estudia los métodosnuméricos de resolución de problemas.

¾Y qué son los métodos numéricos? Un método numérico de resolución de un problemaes un conjunto de reglas que permiten obtener una solución de ese problema (por lo generalno exacta, sino aproximada) mediante la realización de un número �nito de operacioneslógicas y algebraicas elementales. La noción de método numérico se opone por lo general ala de método directo de resolución de un problema, que es aquel que a priori proporcionasoluciones exactas del problema.

Un algoritmo es cualquier descripción ordenada y bien estructurada, paso por paso, de lasoperaciones que conlleva un método �directo o numérico� de resolución de un problema,desde los datos de entrada hasta los resultados. En el ámbito del Cálculo Numérico seestudian también algunos métodos directos, pero en esta asignatura solo trataremos conmétodos numéricos.

Los métodos numéricos requieren normalmente la realización de un número relativamentegrande de operaciones, por lo que el ordenador es una herramienta imprescindible para eldesarrollo y la aplicación del Cálculo Numérico. Los algoritmos que se estudian en estaasignatura se implementarán en el ordenador a través del paquete informático Mathematica.

1Estos ejercicios se indicarán escribiendo entre corchetes la palabra Mathematica al principio del ejercicio.2Hay logros en esta materia desde hace muchos siglos, pero su gran desarrollo comienza a mediados del

siglo XX

5

El coste operacional de un método directo o numérico de resolución de un problema esel número de operaciones que requiere su aplicación. El coste operacional de un método hayque tenerlo en cuenta en relación a la velocidad de cálculo del ordenador y a los errores deredondeo cometidos por el ordenador al realizar las operaciones que conlleva el método. Enel capítulo 1 se precisará la noción de error de redondeo.

Los métodos numéricos son imprescindibles para resolver problemas para los que no seconocen métodos directos de resolución; o bien para problemas para los que existen métodosdirectos pero no son viables (por un excesivo coste operacional, por la importancia de loserrores de redondeo, etc.).

El Cálculo Numérico, ¾se reduce a simples cálculos (aunque sean muchos)?; ¾es una ramamenor de las Matemáticas? La respuesta es en ambos casos negativa: cada método numéricotiene un porqué y unas leyes que se basan en teoremas que tienen tanto peso como encualquier otra rama de las Matemáticas.

Los objetivos generales de la asignatura pueden concretarse en los siguientes.

1. Pretende ser un complemento de otras asignaturas de Matemáticas, introduciendométodos alternativos para la resolución de los problemas que esas asignaturas dejabansin resolver (por ejemplo, la resolución de ciertas ecuaciones o el cálculo de ciertasintegrales que no son alcanzables por métodos clásicos del Álgebra o del Cálculo In-�nitesimal).

2. Se trata de proporcionar una introducción al Cálculo Numérico, mostrando algunos delos métodos numéricos más sencillos y a la vez más utilizados en la resolución de variosproblemas matemáticos importantes.

3. Se procurará que el alumno adquiera una noción de algunas cuestiones importantes quesurgen en el desarrollo y aplicación de los métodos numéricos, tales como el estudiode los distintos errores que pueden producirse, el estudio de la convergencia de losmétodos, etc.

4. Se profundizará en el conocimiento del paquete Mathematica en orden a la imple-mentación de los métodos numéricos estudiados en la asignatura. Se ha elegido elpaquete Mathematica �entre otros paquetes informáticos de cálculo matemático�porque (a) los alumnos ya lo han manejado en otras asignaturas, (b) es un paquetebastante completo y potente, (c) es a la vez relativamente sencillo y (d) está disponibleen las aulas de informática de la Universidad de Almería.

Los objetivos y contenidos especí�cos de la asignatura los detalla su guía docente, quese encuentra en la página web del Departamento de Estadística y Matemática Aplicada dela Universidad de Almería. Estos apuntes explican la teoría de la asignatura y proponenejercicios para desarrollar su parte práctica.

Los conocimientos básicos que se necesitan para cursar con fruto esta asignatura estánincluidos en las asignaturas �Análisis Real�, Álgebra Lineal� y �Matemática discreta� del Plande estudios de la Titulación, si bien solo constituyen una parte menor de esas asignaturasy se repasarán parcialmente en estos apuntes. Lógicamente, también se necesitarán otrosconocimientos matemáticos más básicos pertenecientes a los estudios preuniversitarios.

Agradeceré al lector de estos apuntes que me haga llegar cualquier sugerencia para mejo-rarlos, así como cualquier corrección de erratas o errores que pueda encontrar.

6

Capítulo 1

ANÁLISIS DE ERRORES

En este capítulo se introducen y estudian conceptos y resultados relacionados con dis-tintos tipos de errores numéricos, que suelen producirse en diversos procesos tales como lamedida, los métodos numéricos de que trataremos en los siguientes capítulos y, sobre to-do, como consecuencia de realizar operaciones matemáticas con ordenadores. De hecho, eneste capítulo y en toda la asignatura los ordenadores serán fundamentalmente para nosotrosmáquinas de calcular o computar, calculadoras o computadores que trabajan con números;si bien es cierto que también nos interesará en alguna medida su capacidad para el cálculosimbólico.

1.1. Números máquina

Las calculadoras y los ordenadores realizan las operaciones aritméticas de un modo dis-tinto a como las aprendimos y a como suelen presentarse en los libros de Matemáticas. Paraentender esto, conviene recordar que, si bien los conjuntos de números con los que se sueletrabajar en Matemáticas tienen in�nitos elementos,1 los ordenadores solamente trabajan conun conjunto �nito (aunque pueda ser muy grande) de números, que se denominan números

máquina (cada �máquina� tiene su propio conjunto de números máquina, aunque hoy díaeste conjunto está bastante estandarizado).

En lo que sigue, supondremos que el ordenador con el que trabajamos utiliza el sistema

binario de representación de números; aunque todo lo que se diga se podrá generalizar a otrossistemas de representación. Por brevedad y por simplicidad en la exposición, utilizaremosun modelo un tanto más simpli�cado que los ordenadores más recientes. Esto no supondráperjuicio sino más bien bene�cio para la adquisición de las ideas fundamentales de este capí-tulo. Por otra parte, el estudiante de Informática conoce más a fondo la llamada aritmética

de punto �otante (y su normalización actual) a través de otras asignaturas de su Titulación.

En los ejemplos y ejercicios consideraremos también máquinas que supuestamente traba-jan en el sistema decimal. Esto es útil porque el sistema decimal es más intuitivo, y desdeeste sistema se pueden extrapolar muchas consideraciones al sistema binario; y, a veces, tam-

1Es el caso del conjunto N de los números naturales, del conjunto Z de los números enteros, del conjunto

Q de los números racionales, del conjunto R de los números reales y del conjunto C de los números complejos,

entre otros.

7

bién será útil para resaltar algunas diferencias que existen entre el sistema binario y otrossistemas de representación. Por ejemplo, del mismo modo que en el sistema decimal se tieneque

−2032,5403 = −(2 · 103 + 0 · 102 + 3 · 101 + 2 · 100 + 5 · 10−1 + 4 · 10−2 + 0 · 10−3 + 3 · 10−4),

en el sistema binario se tiene que

−101,11012 = −(1 · 22 + 0 · 21 + 1 · 20 + 1 · 2−1 + 1 · 2−2 + 0 · 2−3 + 1 · 2−4) (= −5,812510).

(a menudo se indica con un subíndice la base o sistema en que está representado un número).Recordemos que tanto en el sistema binario como en el decimal, los números máquina

�excepto el cero� pueden representarse mediante su signo, su mantisa y su exponente. Lamantisa y el exponente vienen determinados por el hecho de que la mantisa, en cualquier base,debe ser de la forma a1.a2a3 . . . de manera que el primer dígito está situado inmediatamentea la izquierda del punto decimal y es no nulo: a1 = 0 (evidentemente, en el caso binarioocurrirá siempre que a1 = 1). Por ejemplo, los dos números anteriores pueden expresarsecomo sigue:

−2032,5403 = −2,0325403 · 103, −101,1101 = −1,011101 · 22.

Así, el signo en ambos casos es `−', la mantisa del primero es 2,0325403 y su exponentees 3, la mantisa del segundo es 1,011101 y su exponente es 2: aunque la mantisa esté enbinario, expresaremos los exponentes �y su base� en base 10; el ordenador, por supuesto,los almacena utilizando el mismo sistema que para la mantisa.

Un ordenador almacena los números utilizando una cadena de bits dividida en tres partes:la primera parte almacena el signo, la segunda almacena el exponente y la tercera la mantisa.La capacidad para representar números de ese ordenador depende, por tanto, de la longitudestablecida para esa cadena de bits, y de la distribución de esos bits entre la mantisa y elexponente (el signo solo necesita un bit). Como consecuencia de esa distribución de los bits,cada ordenador tendrá un rango de valores para la mantisa y un rango de valores para elexponente.

En este capítulo, salvo que se diga otra cosa, supondremos que nuestro ordenador puedeguardar mantisas de hasta 53 dígitos (en el sistema binario); y que el rango para el exponenter de los números máquina sea −210 = −1024 ≤ r ≤ 210−1 = 1023. En los siguientes párrafosvamos a estudiar cuántos números máquina hay en este caso y cómo se distribuyen en larecta real.

En primer lugar, un simple cálculo permite establecer que el conjunto de números máqui-na de ese ordenador �solo� tiene 264+1 elementos (hemos contabilizado al cero como númeromáquina); es decir, hay algo más de 1,8 · 1019 números máquina; muy pocos en relación conlos números reales, que son in�nitos.

En cuanto a cómo se distribuyen los números máquina en la recta real, basta que nos�jemos en la parte positiva de la recta (los números máquina negativos se obtienen cambiandoel signo a los positivos). Si el exponente del número máquina es r = 0, entonces el númerose reduce a su mantisa. Observe que una mantisa de 53 dígitos puede representar hasta252 ≃ 4,5 · 1015 números distintos, que se sitúan en el intervalo [1, 2): el menor de esos

8

números es el 1,02 = 110 y el mayor es 1,11 . . . 112 (en este número hay 53 dígitos iguales a1), que es un valor muy próximo a 210.2 Si el exponente del número máquina es r = 1, denuevo tenemos 252 números máquina, esta vez situados en el intervalo [2, 4); si r = 2, se tieneotros 252 números máquina, todos en el intervalo [4, 8). Por el otro lado, si r = −1 tenemosque hay 252 números máquina en el intervalo [1/2, 1), etc. En general, para todo exponenteentero r tal que −1024 ≤ r ≤ 1023, resulta que hay exactamente 252 números máquina en elintervalo [2r, 2r+1), cuya longitud es 2r. Observe que la densidad de números máquina en lasemirrecta real positiva aumenta a medida que nos acercamos a cero, y disminuye a medidaque nos desplazamos hacia la derecha en dicha semirrecta.

Por tanto, todos los números máquina positivos se encuentran en la unión de todos losintervalos [2r, 2r+1), con −1024 ≤ r ≤ 1023, unión que es igual al intervalo [2−1024, 21024)

(2−1024 ≃ 5,56 · 10−309 y 21024 ≃ 1,80 · 10308). Como consecuencia, el menor número máquinapositivo del ordenador que estamos considerando es precisamente 2−1024 y el mayor es1,11 . . . 112 · 21023 = (2 − 2−52) · 21023 = 21024 − 2971 = 2971(253 − 1) (que también es aproxi-madamente igual a 1,80 · 10308: aunque 2971 sea un número �grande�, resta poco a 21024, quees mucho más �grande�).

Observemos �nalmente que la distancia entre dos números máquina consecutivos dependede la diferencia entre dos valores consecutivos de la mantisa �que en nuestro ordenador esigual a 2−52� y del exponente. Por ejemplo, encontremos los primeros números máquina quesiguen al número 1, que es el número máquina de mantisa 1,0 y exponente 0: 110 = 1,02 · 20.Luego el número máquina siguiente al 1 es el x1 = 1,000 . . . 01 (hay 51 ceros consecutivosen la mantisa), que en base 10 es x1 = 1 + 2−52 ≃ 1 + 2,22 · 10−16.3 Y los números máquinaque siguen al x1 son 1 + 2−52 + 2−52 = 1 + 2−51, 1 + 2−51 + 2−52. . . Por el otro lado, es claroque el número máquina anterior al 1 es x−1 = 1,11 . . . 112 · 2−1 (en la mantisa hay 53 dígitosiguales a 1); es decir x−1 = (2 − 2−52) · 2−1 = 1 − 2−53, que dista de 1 la mitad que x1. Aligual que en estos ejemplos, es fácil concluir que la diferencia entre dos números máquinaconsecutivos situados en el intervalo [2r, 2r+1], con −1024 ≤ r ≤ 1023, es de 2r−52.

En resumen, en el ordenador típico que hemos establecido se producen los siguienteshechos:

1. Hay un total de 264 + 1 ≃ 1,8 · 1019 números máquina (252 ≃ 4,5 · 1015 mantisas,211 = 2048 exponentes y 2 signos, más el número 0).

2. La mantisa menor es 1,00 . . . 0 = 1, y la mayor es 1,11 . . . 1 = 2− 252.

3. Hay un total de 252 números máquina en cada uno de los intervalos [2r, 2r+1) y(−2r+1,−2r] para todo entero r ∈ [−1024, 1023].

4. El menor y mayor números máquina positivos son 2−1024 ≃ 5,56 ·10−309 y 21024−2971 ≃1,80 · 10308, respectivamente; mientras que los negativos son −21024 + 2971 y −2−1024.

5. La distancia entre dos números máquina consecutivos del intervalo [2r, 2r+1] (o delintervalo [−2r+1,−2r]) es de 2r−52, para todo −1024 ≤ r ≤ 1023.

2El número 1,11 . . . 112 es igual al siguiente número en base 10: 20 + 2−1 + 2−2 + · · · + 2−52 = 2 −2−52 = 1,9999999999999997779553950749686919152736663818359375 (número decimal exacto con 52 cifras

decimales, las 15 primeras iguales a 9).3En concreto, x1 es exactamente igual a 1,0000000000000002220446049250313080847263336181640625

(número decimal exacto con 52 cifras decimales, las 15 primeras iguales a 0).

9

Examinando lo expuesto en los párrafos anteriores acerca de la distribución de los númerosmáquina positivos, se pueden observar varias paradojas, como por ejemplo las siguientes:

En el intervalo [0, 2−1024) hay un solo número máquina, el cero; pero en el intervalo[2−1024, 2−1023), contiguo al anterior y de la misma amplitud, hay 252 ≃ 4,5·1015 númerosmáquina.

En el intervalo [2−1024, 2−1023), como hemos dicho, hay 252 números máquina, y es unintervalo de amplitud 2−1024 ≃ 5,56·10−309 (½pequeñísima!). Sin embargo, en el intervalo[21024,∞), que es de amplitud in�nita, no se encuentra ningún número máquina.

Hay tantos números máquina en el intervalo [21023, 21024), que es un intervalo de am-plitud 21023 ≃ 8,99 · 10307 (½enorme!) como en el intervalo [2−1024, 2−1023).

La distancia entre dos números máquina consecutivos en el intervalo [21023, 21024] esde 2971 ≃ 1,996 · 10292; en el otro extremo, la distancia entre dos números máquinaconsecutivos en el intervalo [2−1024, 2−1023) es de 2−1076 ≃ 1,235 · 10−324.

Hay tantos números máquina en el intervalo (0, 1) como en el intervalo [1,∞).

Más adelante volveremos a estas paradojas: en la sección 1.2 veremos cómo subsanaralgunos de los problemas que llevan consigo; y en la sección 1.3 quizá ya no las veamos tanparadójicas.

Con programas como Mathematica es posible trabajar de un modo normal, es decir, conla precisión que permite el uso de los números máquina que acabamos de describir (y asítrabaja el programa, por defecto), pero también es posible aumentar la precisión: veremosesto en las prácticas de la asignatura.

Mencionemos �nalmente que Mathematica sabe trabajar también con números exactoscomo si fueran símbolos, de manera que las operaciones que realiza con ellos son exactas.En este caso, cuando son muchas las operaciones a realizar, lo más probable es que losresultados exactos sean inmanejables y sea necesario obtenerlos de modo numérico (que eslo que ocurrirá habitualmente en este asignatura).

1.2. Redondeo, desbordamiento y cancelación

Una vez que hemos ilustrado cómo son y cómo se distribuyen los números máquina paraun ordenador típico, podremos entender mejor los errores y problemas derivados del hechode que el ordenador solo disponga de esos números para representar cualquier cantidad real.Al introducir números en el ordenador y al operar con ellos, muchas veces tales números o losresultados de las operaciones no son números máquina. Nos preguntamos ahora qué �hace�un ordenador con un número x que no es un número máquina. Pueden darse tres situaciones(las explicamos con respecto a dicho ordenador típico):

1. Que x pertenezca al intervalo [2−1024, 21024 − 2971], donde se encuentran todos losnúmeros máquina positivos; o que x pertenezca al intervalo [2971 − 21024,−2−1024],donde se encuentran todos los números máquina negativos.

2. Que x = 0 pertenezca al intervalo (−2−1024, 2−1024).

10

3. Que x pertenezca al intervalo (21024 − 2971,∞) o al intervalo (−∞, 2971 − 21024).

La respuesta a cada una de las tres situaciones anteriores se da en las tres siguientes subsec-ciones, respectivamente.

1.2.1. Errores de redondeo

Cualquier número que se encuentre en el rango alcanzado por los números máquinapositivos o negativos será redondeado por el ordenador, es decir, sustituido por un númeromáquina. A continuación veremos cómo lo suele hacer.4

El modo usual de redondeo empleado por los ordenadores y calculadoras actuales es elllamado redondeo simétrico, que consiste en lo siguiente. Si x ∈ R no es un número máquinay x ∈ (m1,m2), donde m1 y m2 son números máquina consecutivos, entonces:

1. x se redondea por m1 si x ∈ (m1,m1+m2

2);

2. x se redondea por m2 si x ∈ [m1+m2

2,m2).

Otro sistema de redondeo es el redondeo por corte, que redondea mediante el númeromáquina inferior más próximo: en el caso anterior, cualquier x ∈ (m1,m2) se redondea porm1.

Para cualquier tipo de redondeo, si x∗ es el resultado de redondear un número real x,entonces se de�ne el error de redondeo ex como el valor absoluto de la diferencia entre elnúmero y su redondeo:

ex = |x− x∗|.

En la sección 1.3 extenderemos esta de�nición al redondeo de números complejos, puntosde Rn, etc. El concepto de error de redondeo puede extenderse también al que es producidopor una acumulación de errores de redondeo en el resultado de varias operaciones (véase elejemplo 1.5 más abajo).

Donde se produzcan errores de redondeo normalmente será necesario asegurarse de queesos errores no sobrepasan ciertas cotas, para que el resultado redondeado sea admisible.Esto es común a todo tipo de errores: a lo largo de la asignatura nos ocuparemos de acotarerrores de distinto tipo.

Los dos ejemplos siguientes ilustran los conceptos introducidos sobre el redondeo.

Ejemplo 1.1 Supongamos que una calculadora trabaja en base 10 con ocho dígitos demantisa, y que realiza redondeo simétrico. Como π = 3,14159265358979 . . . , el valor delnúmero π que guarda esta calculadora es π∗ = 3,1415927. Luego el error de redondeo paradicho número es eπ = |π − π∗| ≃ 4,64102 · 10−8. Pero si no tenemos otra informaciónacerca del número π que la que proporciona la calculadora, lo que podríamos deducir acercadel verdadero valor de π sería que se encuentra en el intervalo [3,14159265, 3,14159275):como utiliza redondeo simétrico, el primer valor real cuyo redondeo es π∗ es, efectivamente,a = 3,14159265 y el primer número mayor que a cuyo redondeo ya no es π∗ es b = 3,14159275.Por tanto, eπ = |π − π∗| ≤ |a− π∗| = |b− π∗| = 5 · 10−8. 2

4Los redondeos de números complejos, vectores, matrices, etc. se reducen a los redondeos de los números

reales que los forman. Por ejemplo, para redondear el número complejo a + b i o el vector (a, b), lo que se

hace es redondear cada uno de los números reales a y b.

11

Ejemplo 1.2 Se ha mostrado más arriba que 1 = 1,00 . . . 00 · 20 y 1 + 2−52 = 1,00 . . . 01 · 20(hay 53 dígitos en las mantisas) son dos números máquina consecutivos en nuestro ordenadortipo. El punto medio entre esos números es 1+1+2−52

2= 1 + 2−53 = 1,00 . . . 001 · 20 (hay 54

dígitos en esta mantisa: no es un número máquina).

Si el ordenador utilizara el redondeo por corte, entonces todos los números x del intervalo(1, 1+2−52) los redondea por x∗ = 1, y el error de redondeo satisface que ex = |x−x∗| < 2−52.En general, cualquier número cuya mantisa binaria tenga más de 53 dígitos, el redondeo porcorte se queda con los primeros 53 dígitos y elimina los demás.

Si el ordenador utiliza el redondeo simétrico, que es lo habitual, entonces todo númerox del intervalo (1, 1 + 2−53) lo redondea por x∗ = 1; mientras que cualquier número x delintervalo [1 + 2−53, 1 + 2−52) lo redondea por x∗ = 1 + 2−52. En uno y otro caso es claroque el error de redondeo puede acotarse como sigue: ex = |x− x∗| ≤ 2−53. En general, paracualquier número cuya mantisa binaria tenga más de 53 dígitos, si el 54o dígito es igual a0, entonces el redondeo simétrico �como el redondeo por corte� se queda con los primeros53 dígitos y elimina los demás; y si el 54odígito es igual a 1, entonces el redondeo simétricoelige el número máquina siguiente al que proporciona el redondeo por corte. 2

Para cualquier ordenador o calculadora, el épsilon de la máquina se de�ne como el menornúmero positivo ε tal que la suma numérica 1 + ε da como resultado un número (númeromáquina) mayor que 1. Por tanto, si el ordenador del ejemplo 1.2 utiliza el redondeo simétrico,entonces el épsilon de la máquina es ε = 2−53 ≃ 1,11 · 10−16, ya que cuando el ordenadorhace la suma 1 + η, para cualquier número η ∈ (0, 2−53), da como resultado 1, mientrasque la suma 1 + 2−53 da como resultado 1 + 2−52 > 1. Si el ordenador aplicara el redondeopor corte, entonces el épsilon de la máquina sería ε = 2−52 (el doble que para el redondeosimétrico): para cualquier η ∈ (0, 2−52), el ordenador suma 1 + η = 1; y la suma 1 + 2−52 dacomo resultado 1 + 2−52 > 1.

Para un ordenador o calculadora dados, es muy conveniente conocer qué precisión puedeesperarse de los redondeos que realiza. La cuestión más simple en este sentido es saber cuántosdígitos correctos consecutivos tendrá la mantisa (en base 10) del redondeo de un númeroreal cualquiera, contados a partir del primer dígito. Por ejemplo, en el ordenador modelo quevenimos usando, habrá unos 15 o 16 dígitos correctos, tanto con redondeo simétrico comocon redondeo por corte. Veamos un par de ejemplos.

Ejemplo 1.3 Partimos del ejemplo 1.2. Si y = 1 + 2−53 = 1,00000000000000011102 . . . , suredondeo simétrico es y∗s = 1 + 2−52 = 1,00000000000000022204 . . . y su redondeo por cortees y∗c = 1. Por tanto, los redondeos y∗s e y∗c coinciden con y en los primeros 16 dígitos de sumantisa, que son: 1,000000000000000. 2

Ejemplo 1.4 En el ordenador habitual, hemos mostrado más arriba que el número máquinaanterior al 1 es 1−2−53. Consideramos el punto medio z entre 1−2−53 y 1, que no es un númeromáquina. Es claro que z = 1− 2−54 = 0,99999999999999994448 . . . , y que sus redondeos porcorte y simétrico son, respectivamente, z∗c = 1 − 2−53 = 0,9999999999999998889 y z∗s = 1.Es claro que z∗c coincide con z en los primeros 15 dígitos de su mantisa (en base 10), queson: 9,99999999999999; teniendo en cuenta que 1 = 0,999 . . . (in�nitos dígitos iguales a 9),

12

y tomando 9,999 . . . como mantisa del número 1, entonces z∗c coincide con z en los primeros16 dígitos de su mantisa, que son: 9,999999999999999. 2

El lector interesado puede encontrar nociones más rigurosas del concepto de precisiónen la bibliografía, estudiando por ejemplo el concepto de dígitos signi�cativos (véanse porejemplo [2] y [11]), que se basa en el concepto de error relativo, que estudiaremos en la sección1.3. Por otra parte, veremos que el error relativo producido por el redondeo es también unamedida de la precisión del redondeo.

Los errores de redondeo se producen tanto en la entrada y en la salida de los datos comoen todas las operaciones intermedias, y pueden acumularse o cancelarse. La limitación quepara los ordenadores supone la producción de errores de redondeo es aceptable la mayorparte de las veces, pero en ocasiones puede ocurrir que los redondeos des�guren de maneraimportante el resultado �nal de las operaciones. En esta asignatura aprenderemos algunosmodos concretos de eludir este problema en determinadas situaciones. Sin embargo, el estudiogeneral de la propagación de los errores de redondeo está fuera de los objetivos de este curso.Una introducción se puede encontrar en muchos libros de Análisis Numérico: por ejemplo,una muy sencilla se encuentra en [9].

Para terminar este apartado, consideramos un ejemplo en el que se realizan un par deoperaciones con el programa Mathematica y se estudia el error de redondeo que se produceal realizarlas.

Ejemplo 1.5 En Mathematica,√3. representa la raíz cuadrada numérica de 3. Por tanto,

x∗ =√3. es un redondeo de x =

√3 (la raíz cuadra exacta). Con la precisión y el redondeo

que utiliza por defecto, Mathematica produce el siguiente resultado:√3. = 1,7320508075688772.

Si Mathematica presenta el resultado redondeado simétricamente a 6 dígitos, como hacepor defecto, entonces muestra

√3. = 1,73205, pero el valor que guarda de

√3. es el que

hemos expresado; o mejor dicho, el número en binario cuya representación en base 10 esaproximadamente igual a 1,7320508075688772.

Teniendo en cuenta que el valor exacto de√3 es

√3 = 1,732050807568877293527 . . . ,

resulta que el error de redondeo es e√3 = |√3−

√3.| = 9,3527 . . . · 10−17 y las mantisas (en

base 10) de√3 y

√3. coinciden en sus primeros 17 dígitos.

Si realizamos una segunda operación, elevando al cuadrado√3., Mathematica obtiene

(√3.)2 = 2,9999999999999996

(hay 15 dígitos iguales a 9). El resultado exacto de (√3.)2 es 3, por lo que el error de

redondeo que se ha producido al realizar las dos operaciones es e3 = |3− (√3.)2| ≃ 4 ·10−16.5

Tomando 2,999 . . . (con in�nitos dígitos iguales a 9) como mantisa del resultado exacto 3,resulta que dicha mantisa coincide con la del redondeo (

√3.)2 en sus primeros 16 dígitos,

que son 2,999999999999999. 2

5Hemos escrito �≃ 4 · 10−16� en vez �= 4 · 10−16� porque, como se ha dicho, el número en base 10

que muestra Mathematica no suele ser exactamente igual al correspondiente número en base 2 que guarda el

programa. De hecho, la respuesta que da Mathematica a la operación 3−(√3.)2 es 4,440892098500626 ·10−16

y no 4 · 10−16.

13

1.2.2. Errores de cancelación

Estudiamos ahora la segunda situación, que se produce cuando el ordenador recibe oproduce un número real x no nulo que pertenece al intervalo (−2−1024, 2−1024) (nos referimosal ordenador habitual). En este caso se produce una situación que suele denominarse deunder�ow : el ordenador redondea el número x por el número 0.

Un problema mucho más grave se produce cuando, a resultas de acumulación de erroresde redondeo, cierto resultado, quizá próximo a 0 (pero no tanto como para situarse en elintervalo (−2−1024, 2−1024): recuerde que 2−1024 ≃ 5,56 · 10−309), o incluso ni siquiera próximoa 0 (de modo absoluto; de modo relativo sí será proximo a 0), se redondea por el número0. Entonces se dice que se ha producido una cancelación. Estas cancelaciones pueden sermuy peligrosas, sobre todo cuando se dan en operaciones intermedias de un proceso: si nose descubren pueden producir errores muy graves en los resultados �nales.

Para uni�car términos, consideraremos que un under�ow es un tipo particular de can-celación (el problema en sí es el mismo: una cantidad no nula que se redondea por 0), de ahíque ese sea el título dado a esta subsección.

Un caso frecuente de cancelación surge cuando se calcula la diferencia de dos númerosbastante próximos, y el redondeo del ordenador iguala esa diferencia a 0. Por ejemplo, six = 1+2−54 e y = 1, nuestro ordenador establecerá que x−y = 0 (ya que, si usa la precisiónhabitual, x = 1+2−54 = 1 = y), cuando en realidad x− y ≃ 5,5 · 10−17. Otro ejemplo similarse obtiene sumando un número natural a los exponentes (1 = 20): si u = 260 + 26 y v = 260,entonces el ordenador calculará u− v = 0, cuando en realidad u− v = 64.

Para evitar los errores de cancelación se pueden emplear diversas tácticas según el tipode cancelación. Son tácticas que son útiles para ciertas situaciones particulares. Aquí soloexaminaremos dos de ellas, ambas aplicables a la evaluación de funciones en ciertos valoresen los que se produce una cancelación: es el caso de una función f tal que f(x) = 0 paracierto valor x de su dominio, pero el ordenador, al redondear, produce como resultado quef(x) = 0.

La primera táctica se aplica cuando f realmente se anula en algún punto a relativamentepróximo al punto x donde f se cancela. En este caso, quizá pueda evitarse la cancelaciónhaciendo el desarrollo de Taylor de la función a evaluar f respecto del punto a donde seanula. Recordamos a continuación el teorema que proporciona el mencionado desarrollo deTaylor. En su enunciado se utiliza la siguiente de�nición.

Sea k ∈ N y sea f una función real de�nida en un intervalo I. Se dice que f es de claseCk en I si f es k veces derivable en I y además su derivada k-ésima f (k) es continua en I.

El teorema mencionado es el siguiente.

Teorema 1.1 (de Taylor) Sea k un número natural, y sea f una función de clase Ck en

un intervalo I tal que existe f (k+1)(x) para todo punto x del interior de I. Entonces, para

cualesquiera a, x ∈ I, con x = a, existe un punto cx entre a y x tal que

f(x) = f(a) +k∑

i=1

f (i)(a)

i!(x− a)i +

f (k+1)(cx)

(k + 1)!(x− a)k+1. (1.1)

A la expresión (1.1) se le denomina fórmula o desarrollo de Taylor de orden k de lafunción f en el punto a.

14

En los problemas propuestos mostramos cómo se aplica este teorema para evitar cancela-ciones.

Una segunda táctica que suele emplearse cuando aparece una diferencia de funciones quese cancelan es multiplicar y dividir la expresión por su conjugado (por la suma de esas dosfunciones). Veamos un ejemplo de esta situación.

Ejemplo 1.6 Al utilizar una calculadora pequeña o incluso Mathematica, tendremos unproblema al calcular valores de la función

f(x) =1√

x+ 1−√x

(1.2)

en valores grandes de la variable x. Por ejemplo, para valores de la forma x = 10i, con i unnúmero natural positivo su�cientemente grande (la cancelación produce en este caso un errorde la forma 1/0). Sin embargo, multiplicando en f(x) por el conjugado del denominador seevita la cancelación, ya que se obtiene que

f(x) =√x+ 1 +

√x.

Pruebe lo que pasa en su calculadora o con Mathematica al calcular por ambos métodosf(1010), f(1015) o f(1020). Por ejemplo, al tomar x = 1020 en (1.2) se producirá un error, yaque el denominador será redondeado por cero. Sin embargo, al usar la segunda expresión seobtendrá que f(1020) =

√1020 + 1 +

√1020 ≃ 2 ·

√1020 = 2 · 1010. 2

1.2.3. Errores de desbordamiento

Estudiamos ahora la tercera situación considerada al principio de esta sección, que seproduce cuando el ordenador (nos referimos al ordenador típico) recibe o produce un númeroreal x que pertenece al intervalo (21024 − 2971,∞) o al intervalo (−∞, 2971 − 21024). En estecaso se produce una situación que suele denominarse de over�ow o desbordamiento, que elordenador indicaría de alguna manera (con las palabras over�ow, error, etc.).

Para evitar los errores producidos por desbordamiento se pueden emplear diversas tác-ticas, que dependen mucho del caso concreto. Como ocurre con la cancelación, el desbor-damiento suele producirse en los pasos intermedios de una serie de operaciones, mientras queel resultado �nal de ellas sea perfectamente representable por el ordenador. En ese caso, sesuele buscar un reescalamiento (un cambio de escala) de las operaciones (por ejemplo, sacan-do factor común una de las variables que contribuye a que se produzca el desbordamiento).

Veamos un ejemplo ilustrativo, esta vez �por simplicidad� propio de una calculadorapequeña (usando Mathematica no habría problemas de desbordamiento para este ejemplo).

Ejemplo 1.7 Si tiene una calculadora pequeña, es muy posible que no pueda calcular di-rectamente la diferencia 33144 − 101500, ya que 33144 y 101500 son números demasiado grandespara la calculadora. Pero sí puede aproximar bien el resultado si realiza el cálculo de lasiguiente forma:

33144 − 101500 = 101500(

33144

101500− 1

)= 101500

([33,144

101,5

]1000− 1

)≃ 101500

(1,0001594087861000 − 1

)≃ 0,17280238 · 101500.

15

Por supuesto, el resultado se sitúa fuera del rango de la calculadora (e incluso del de nuestroordenador típico), pero es manejable y se tiene una idea clara de su magnitud. 2

1.3. Error absoluto y error relativo

En numerosas ocasiones, el error (de redondeo, de medida, de aproximación, etc.) quepuede presentar un número, un punto, un vector, etc. suele describirse desde uno de los puntosde vista que ofrecen los dos conceptos que se estudian en esta sección y que introducimos acontinuación.

De�nición 1.1 Sea x un número real y sea x∗ una aproximación de ese número. Entonces,el error absoluto ex de la aproximación se de�ne como el siguiente número:

ex = |x− x∗|; (1.3)

y �siempre que x = 0� su error relativo rx es la siguiente proporción:

rx =ex|x|

=

∣∣∣∣x− x∗

x

∣∣∣∣ . (1.4)

Por ejemplo, el error de redondeo de�nido en la subsección 1.2.1 es el error absoluto dela aproximación de un número por su redondeo. Veamos un ejemplo en el que también secalcula el error de redondeo relativo.

Ejemplo 1.8 En la sección 1.1 se ha mostrado que �para el ordenador típico establecido�la distancia entre dos números máquina consecutivos del intervalo [2r, 2r+1] es de 2r−52, paratodo −1024 ≤ r ≤ 1023. Y que con el redondeo simétrico, el punto medio de esos intervalosse redondea por el número máquina mayor. Así, en el ejemplo 1.2, para r = 0, mostramosque si x = 1 + 2−53 su redondeo era x∗ = 1 + 2−52. Si consideramos r = 60 y u = 260 + 27,resulta que u ∈ [260, 260 + 28], que es un intervalo cuyos extremos son dos números máquinaconsecutivos. Es claro que el redondeo de u es u∗ = 260+28. Por tanto, los errores de redondeopropiamente dichos (es decir, los absolutos) son

ex = |x− x∗| = 2−53 ≃ 1,11 · 10−16 eu = |u− u∗| = 27 = 128.

Sin embargo los errores relativos son

rx =ex|x|

=2−53

1 + 2−53ru =

ex|x|

=27

260 + 27=

2−53

1 + 2−53.

Luego la distribución de números máquina a lo largo de la recta real ya no es tan paradóji-ca, puesto que genera errores relativos similares. 2

En el ejemplo anterior se observa el interés del error relativo. Este informa mejor queel absoluto acerca de la precisión de una aproximación, ya que es una proporción de errorrespecto del número aproximado. Evidentemente, no es lo mismo que haya un error absolutode 2 kg en una báscula para personas que en una para camiones.

16

El error relativo suele darse también en tanto por ciento: un error relativo de rx (que esun �tanto por uno�, una proporción, como hemos dicho) es lo mismo que un error relativodel 100 · rx%.

Hay libros que de�nen los conceptos de la de�nición 1.1 sin valor absoluto, pero de hecholuego utilizan principalmente el valor absoluto de los errores. Además, con valor absoluto sefacilita la extensión de dichos conceptos a otros conjuntos: la de�nición 1.1 puede extendersea Rn o, en general, a cualquier espacio normado, concepto que introducimos a continuación.

De�nición 1.2 Un espacio normado V es un espacio vectorial real en el que se ha de�nidouna norma, es decir, una aplicación ∥ · ∥ : V → [0,∞) tal que, para cualesquiera u, v ∈ V yα ∈ R, se satisfacen las siguientes propiedades:

(a) ∥u∥ = 0 si y solo si u es el elemento neutro de V ;

(b) ∥αu∥ = |α|∥u∥;

(c) ∥u+ v∥ ≤ ∥u∥+ ∥v∥.

El espacio euclídeo n-dimensional es un ejemplo de espacio normado. Recuerde que lanorma euclídea de un vector u = (u1, u2, . . . , un) de Rn se de�ne como sigue:

∥u∥ =√

u21 + u2

2 + · · ·+ u2n.

Observe que las tres propiedades de la de�nición de norma tienen una sencilla interpretacióngeométrica en el caso de la norma euclídea en R2 y R3.

Otro ejemplo de espacio normado es el conjunto C de los números complejos, tomandocomo norma el módulo.

En todo espacio normado V se puede de�nir una distancia o métrica d : V ×V → [0,∞),de�nida por

d(u, v) = ∥u− v∥ para cualesquiera u, v ∈ V.

En un espacio normado (incluido R), las fórmulas (1.3) y (1.4) de los errores absoluto yrelativo se escriben respectivamente de la siguiente forma:

ex = ∥x− x∗∥, rx =ex∥x∥

=∥x− x∗∥

∥x∥.

Observe que el error absoluto es la distancia entre x y su aproximación x∗: ex = d(x, x∗).En muchos problemas del Cálculo Numérico interesa tener alguna información acerca de

la magnitud del error absoluto o del error relativo que se ha producido en un resultado.Ya lo hemos visto en el caso del redondeo. Como en este caso, dicha información consistirámuchas veces en la obtención de cotas de esos errores. Lógicamente, se procurará que esascotas sean tan ajustadas como se pueda conseguir: cuanto menores sean, mejor informaciónproporcionan sobre los errores. Tanto en este capítulo como en los siguientes, cuando se pidaen algún ejercicio que se obtenga una cota de un determinado error, se sobreentenderá, portanto, que dicha cota debe ser la menor que se pueda conseguir. Del mismo modo, cuandose pida un intervalo donde se encuentre un determinado número que se está aproximando(podrá ser una integral, la solución de una ecuación, etc.), también se entenderá que debe

17

obtenerse a partir de las mejores cotas que se puedan encontrar. E igualmente se actuará enotros problemas similares.

En el ejemplo 1.1 obtuvimos un cota de un error de redondeo (absoluto). En concreto,probamos que eπ ≤ 5 · 10−8. En cuanto al error relativo, se tiene que

rπ =eππ

≤ 5 · 10−8

π<

5 · 10−8

3,14159265=

5

314159265≃ 1,59155 · 10−8.

El siguiente resultado generaliza lo que acabamos de realizar para acotar rπ, puesto quepermite acotar el error relativo cuando se tiene una cota del error absoluto, siempre que sesatisfaga una hipótesis bastante liviana.

Teorema 1.2 Supongamos que x∗ es una aproximación de un elemento x de un espacio

normado (V, ∥ · ∥). Si se conoce una cota εx del error absoluto ex tal que εx < ∥x∗∥, entoncesse tiene que

rx ≤ εx∥x∗∥ − εx

.

Observe que la condición εx < ∥x∗∥ es casi irrelevante: es un caso muy extraño que lacota del error absoluto sea mayor o igual que el valor absoluto de la aproximación.

El siguiente resultado permite acotar el error absoluto cuando se tiene una cota del errorrelativo, siempre que esta cota cumpla una condición muy poco exigente.

Teorema 1.3 Supongamos que x∗ es una aproximación de un elemento x de un espacio

normado (V, ∥ · ∥). Si se conoce una cota ρx < 1 del error relativo rx, entonces se tiene que

ex ≤ ρx1− ρx

∥x∗∥.

1.4. Ejercicios

1. Supongamos que se dispone de un procesador muy simple, que trabaja en base diez,con una cifra de mantisa y exponente entre −1 y 1.

a) Obtenga todos los números máquina que puede manejar este procesador.

b) Represente los números máquina anteriores en una recta real. Entre ellos apare-cerán huecos: ¾son todos los huecos de la misma longitud?

c) Obtenga los errores absolutos y relativos que comete este procesador al representarlos números 1,3, 0,13, 1,02 y 0,102. Obtenga los errores relativos de dos formas,una de ellas utilizando el teorema 1.2

2. [Mathematica] Muestre en la pantalla de su ordenador todos los números máquina deun procesador decimal de dos cifras de mantisa y exponente entre −3 y 3.

3. [Mathematica] Encuentre el épsilon de la máquina donde trabaja (para lo cual, supongaque este épsilon será un número de la forma 2−n, con n ∈ N). Compare el épsilon de lamáquina con los menores números de la forma 2−n tales que en Mathematica resulteque (a) 10−6 + 2−n sume más que 10−6 y (b) 106 + 2−n sume más que 106.

18

4. Las calculadoras Casio de la serie fx trabajan en base diez, con 10 cifras de mantisa yexponente entre −98 y 100. Supongamos que usamos una de esas calculadoras.

a) ¾Cuáles son los dos números máquina consecutivos entre los que se halla el númeroe = 2,718281828459 . . .? ¾Cómo redondea la calculadora el número e si hace re-dondeo simétrico? Acote el error de redondeo que se produce en este caso.

b) ¾Cuál es el número máquina que sigue al número 4,38938 · 102? ¾Y cuál es el quesigue a 0,000003289?

c) ¾Existe la misma distancia entre cada uno de los tres números máquina consi-derados y su siguiente? Encuentre dos números máquina consecutivos que distenentre sí exactamente 10−8.

5. Para cada uno de los números reales siguientes, halle el error absoluto y el error relativode su aproximación (que se indica con un superíndice �∗�).

a) x = 2,71828182 x∗ = 2,7182

b) y = 98350 y∗ = 98000

c) z = 0,000068 z∗ = 0,00006

6. Las primeras cifras de√2 son las del número r = 1,4142135623. Si solo se dispone

de este dato acerca de√2 y se aproxima esta raíz cuadrada por el número r, obtenga

una cota del error absoluto que se comete al realizar esa aproximación. Obtenga unaacotación mejor de ese error absoluto |

√2 − 1,4142135623| si se sabe que la siguiente

cifra decimal de√2 es 7. Repita lo mismo si se conociera una cifra más, esto es, que√

2 ≃ 1,414213562373.

7. Se dispone de una calculadora que trabaja en base diez, con 7 cifras de mantisa yexponente entre −25 y 25.

a) Indique entre qué dos números máquina consecutivos se encuentra el número√6

(= 2,449489742783 . . .).

b) Si la calculadora utiliza el redondeo por corte, ¾con qué número máquina identi-�caría al número

√6? ¾Y si utiliza el redondeo simétrico?

c) En el caso del redondeo simétrico, y suponiendo que solo dispone de la calculado-ra del ejercicio, obtenga: (c1) las expresiones de los errores absoluto y relativocometidos; (c2) una cota de dichos errores.

8. Se ha medido la altura de un edi�cio, resultando que esta es de 83,72 metros. Todoslos dígitos de esta cantidad se conocen con certeza, pero no es posible precisar ningunode los siguientes. Obtenga un intervalo de números reales en el que se pueda asegurarque se encuentra la altura exacta del edi�cio. Acote el error relativo que se comete alaproximar la altura exacta por el valor medido.

9. La temperatura de un �uido medida por un termómetro es de 224,3◦K. Se sabe queel termómetro comete como máximo un error relativo del 0,02%. Calcule una cotadel error absoluto cometido en la medida de la temperatura del �uido y determine unintervalo donde se pueda asegurar que se encuentra la medida correcta.

19

10. La solución exacta de un problema es (245, 310, 125). Un método de aproximación dela solución produce como resultado (243, 311, 124). ¾Cuál es el error absoluto y el errorrelativo cometido por dicho método?

11. Considere el cálculo de la función f(x) =√x3 + 4− 2 para x próximo a cero. Muestre

que para valores su�cientemente pequeños de x se produce la cancelación del valorde f(x). Compruebe que esta situación puede evitarse si se multiplica y divide por elconjugado de f(x) para obtener que f(x) = x3/(

√x3 + 4 + 2). Compruebe también

que la cancelación también se evita si se utiliza el desarrollo de Taylor de la función f .

12. a) Encuentre un modo de evaluar las siguientes expresiones en los valores que seindican, de modo que puedan evitarse errores de cancelación.1) 1− cos x en valores de x próximos a cero.2)

√x4 + 4− 2 en valores de x próximos a cero.

3) x(√x+ 1−

√x) en valores grandes de la variable x.

4) x− 1/2−√x2 − x+ 3 en valores grandes de la variable x.

5) 2−√x2 − 6x+ 13 en valores próximos a x = 3.

b) [Mathematica] Por otra parte, para cada uno de los apartados anteriores, utiliceMathematica para construir una tabla en la que se comparen los resultados dela evaluación de (1) la expresión dada y de (2) la expresión obtenida para evitarla cancelación, en una sucesión de números que tienda al valor cerca del cual seproduce esa cancelación.

13. ¾Para qué valores de x se pueden presentar problemas de cancelación en el cálculo deex − e y de ln x− 1? Encuentre un modo de resolver dichas situaciones.

14. Acote el error absoluto y el error relativo que comete Mathematica al aproximar√2

por√2. = 1,4142135623730951.

15. Una calculadora trabaja en base diez, con 12 cifras de mantisa y exponente entre−120 y120, y aplica redondeo simétrico. Consideramos el número x = 718,2288118822884455.

a) ¾Entre qué dos números máquina consecutivos se encuentra x?

b) Obtenga los errores absoluto y relativo que comete la calculadora al aproximar elvalor de x por el número máquina correspondiente.

16. Se dispone de una calculadora que muestra hasta 12 dígitos decimales y utiliza redondeosimétrico. Si al realizar un cálculo muestra −12345678, 9012 como resultado, ¾qué sepuede deducir acerca del valor exacto del cálculo? Proporcione una cota del errorabsoluto y del error relativo que ha cometido la calculadora al mostrar dicho resultado.

17. Una calculadora trabaja en base diez, con 10 cifras de mantisa y exponente entre −100

y 99, y aplica redondeo por corte.

a) El número x = 0,00123456789876, ¾entre qué dos números máquina consecutivosse encuentra?

b) Obtenga los errores absoluto y relativo (expresados en notación cientí�ca: signo,mantisa y exponente) que comete la calculadora al aproximar el valor de x por elnúmero máquina correspondiente.

20

18. El peso medido en una báscula de precisión de una porción de mineral es de 48,125

gramos. El error máximo de la báscula es del 0,2%. Determine un intervalo donde sepueda asegurar que se encuentra el peso exacto de esa porción de mineral.

19. Un mortero dispara repetidas veces contra un blanco situado en una zona llana. Elblanco está a 300 metros del mortero. Se estima que el error máximo que puede cometeren circunstancias normales es de 12 metros alrededor del blanco. Obtenga una cota delerror relativo que puede tener el disparo del mortero.

Algunas cuestiones teóricas

1. De�na los conceptos de número máquina y de épsilon de una máquina.

2. Explique brevemente qué es un error de cancelación y las consecuencias que puedetener.

21

Capítulo 2

RESOLUCIÓN NUMÉRICADE ECUACIONES Y SISTEMAS NOLINEALES

En este capítulo nos interesaremos por los métodos de localización y aproximación de lassoluciones de una ecuación de una incógnita real (o compleja, en algunos casos). Tambiénharemos una breve introducción a la aproximación de soluciones de sistemas de ecuacionesno lineales.

2.1. Introducción a los números complejos

Dedicamos esta primera sección a introducir algunas nociones básicas sobre númeroscomplejos. Hay varias formas de introducir estos números. Aquí nos limitamos a introducirlossin gran rigor matemático, buscando solo su utilidad como herramienta.

Como estamos en un capítulo dedicado a resolver ecuaciones, vamos a introducir losnúmeros complejos como soluciones de un tipo de ecuaciones. Sabemos que las solucionesde la ecuación x2 = 1 son x = 1 y x = −1. Ambos números se representan en la recta realmediante un punto cuya distancia al origen 0 es igual a 1. Sabemos también que la ecuaciónx2 = −1 no tiene soluciones reales (ningún número real al cuadrado es igual a −1). Entoncesse de�ne la unidad imaginaria i como la raíz cuadrada de −1:

i =√−1.

También se puede considerar el número imaginario opuesto − i = −√−1. Operando como

en R resulta que i2 = (− i)2 = −1, de modo que ± i son las dos soluciones (veremos que nopuede tener más) de la ecuación x2 = −1.

Entonces, se de�ne el conjunto C de los números complejos como el formado por losnúmeros de la forma

a+ b i, para cualesquiera a, b ∈ R.

Observe que R ⊂ C, ya que los números reales son los complejos de la forma a + 0i,con a ∈ R. A los números complejos que no son reales �es decir, los complejos a + b i

tales que b = 0� se les llama números imaginarios. A los números complejos de la forma

23

0 + b i = b i, con b ∈ R − {0}, se les llama números imaginarios puros. Observe que losnúmeros imaginarios puros son las soluciones de todas las ecuaciones de la forma x2 = −b2,con b ∈ R − {0}. En coherencia con las de�niciones anteriores, a los números a y b se lesllama parte real y parte imaginaria del número complejo a+ b i, respectivamente.

Los números complejos se pueden representar en el plano xy (que en este contexto sueledenominarse plano complejo), de modo que cada número complejo a + b i se representamediante el punto (a, b) (la parte real del número complejo es la abscisa del punto quelo representa, mientras que la parte imaginaria es la ordenada). En particular, los númerosreales se representan en el eje x, y los imaginarios puros en el eje y.

En los tres párrafos siguientes introducimos tres de�niciones importantes que puedenaplicarse �con una excepción, que mostraremos� a cualquier número complejo z = a+ b i.

El conjugado de z se de�ne como el número complejo z = a − b i. Observe que z y z

son puntos simétricos respecto del eje x en el plano complejo (tienen la misma abscisa yordenadas opuestas).

El módulo de z se de�ne como el número real |z| =√a2 + b2. Observe que el módulo de

a + b i coincide con la distancia del punto (a, b) al punto (0, 0) �que coincide a su vez conel módulo o norma del vector (a, b) de R2� y que, por tanto, |0| = 0 y |z| > 0 si z = 0. Esigualmente trivial que |z| = |z| para todo z ∈ C.

El argumento principal de un número complejo no nulo z = a + b i, que denotamos porarg z, se de�ne como el ángulo que forma el vector de posición del punto (a, b) con el semiejepositivo de abscisas. Luego arg z ∈ [0, 2π). El argumento principal no está de�nido, portanto, para z = 0.

En la �gura 2.1 se representa un número complejo z así como su módulo y su argumentoprincipal en el plano complejo.

z

z

| |z

arg

Figura 2.1: Módulo y argumento principal de un número complejo z

Si z = a + b i y denotamos θ = arg z, entonces es claro que a = |z| cos θ y b = |z| sen θ.Por tanto, z = a + b i = |z| cos θ + |z| sen θ · i = |z|(cos θ + i sen θ). Mientras que a a + b i sele llama forma binómica del número complejo z, a

z = |z|(cos θ + i sen θ), donde θ = arg z,

se le llama forma trigonométrica de z. Esta fórmula es válida para el número z = 0 si seconsidera que su argumento principal es cualquier valor del intervalo [0, 2π).

24

Evidentemente, para cualquier z ∈ C, también es cierto que z = |z|(cos θ + i sen θ) siθ = arg z + 2kπ, con k ∈ Z− {0}. A estos valores de θ se les llama también argumentos dez: por eso se dice que arg z es el argumento principal de z.

El argumento principal de un número complejo z = a+ b i = 0 se puede calcular de unamanera sencilla teniendo en cuenta los signos de a y b y que tg θ = b/a si a = 0, obteniéndoseel siguiente resultado.

1. Si a > 0 y b ≥ 0, entonces arg z = arc tg(b/a).

2. Si a > 0 y b < 0, entonces arg z = arc tg(b/a) + 2π.

3. Si a = 0 y b > 0, entonces arg z = π/2.

4. Si a = 0 y b < 0, entonces arg z = 3π/2.

5. Si a < 0, entonces arg z = arc tg(b/a) + π.

En el conjunto de los números complejos se de�nen las operaciones suma y producto, delas que se derivan la diferencia y el cociente, que extienden las operaciones respectivas enel conjunto de los números reales, así como sus propiedades básicas. Las tres primeras sede�nen de la siguiente forma:

(a+ b i) + (c+ d i) = a+ c+ (b+ d)i

(a+ b i)− (c+ d i) = a− c+ (b− d)i

(a+ b i)(c+ d i) = ac− bd+ (ad+ bc)i.

La fórmula del producto se puede deducir de la propiedad distributiva, lo que evita tenerque memorizarla:

(a+ b i)(c+ d i) = ac+ ad i + bc i + bd i2 = ac− bd+ (ad+ bc)i.

En cuanto al cociente de dos números complejos, es trivial el cociente o división de unnúmero complejo entre un número real no nulo:

a+ b i

c=

a

c+

b

ci.

Pero para obtener el cociente de un número complejo entre un número imaginario se requieremultiplicar numerador y denominador por el conjugado del denominador, con lo que se llegaa un nuevo denominador que es un número real no nulo, y entonces se opera como acabamosde hacer:

a+ b i

c+ d i=

(a+ b i)(c− d i)

(c+ d i)(c− d i)=

ac+ bd+ (−ad+ bc)i

c2 + d2=

ac+ bd

c2 + d2+

bc− ad

c2 + d2i.

Con las operaciones realizadas en el denominador de la fórmula anterior, se ha probadotambién que el producto de todo número complejo por su conjugado coincide con el cuadradode su módulo: zz = |z|2 para todo z ∈ C.

Las potencias enteras de un número complejo se de�nen igual que para los númerosreales, y tienen propiedades semejantes. Un caso interesante es el de las potencias de la

25

unidad imaginaria i. Por de�nición i0 = 1, i1 = i e i2 = −1. Por tanto, i3 = i2 i = − i ei4 = (i2)2 = 1 = i0. Luego

1 = i0 = i4 = i8 = · · · ,i = i1 = i5 = i9 = · · · ,

−1 = i2 = i6 = i10 = · · · ,− i = i3 = i7 = i11 = · · · .

En general, para calcular in, con n ∈ N, se elevará i al resto de la división de n entre 4, ylos únicos resultados de estas potencias son 1, i, −1 y − i. Y para calcular i−n basta escribiri−n = 1/ in y efectuar la división. Los únicos resultados posibles de las potencias negativassiguen siendo los mismos que los de las positivas.

2.2. Soluciones de ecuaciones y su multiplicidad

En esta sección recordaremos la noción de multiplicidad de una raíz de un polinomio (ofunción polinómica1) y la extenderemos a cualquier tipo de función o de ecuación.

La forma general de una ecuación con una incógnita x es g(x) = h(x), donde g y h

son funciones. Pero siempre es posible pasar h(x) al primer miembro o g(x) al segundo,resultando una ecuación de la forma

f(x) = 0. (2.1)

A lo largo de este capítulo la ecuación a resolver la expresaremos casi siempre de esta forma(una función igualada a cero), lo que facilitará enormemente la introducción de los distintosconceptos y resultados que estudiaremos.

El problema de encontrar las soluciones (también se les llama raíces) de una ecuaciónf(x) = 0 es evidentemente equivalente al problema encontrar los ceros (también llamadosraíces) de la función f(x). Observe que el término �raíz� lo hemos usado indistintamentepara hablar de las soluciones de f(x) = 0 o de los ceros de f .

2.2.1. Ecuaciones polinómicas, raíces y su multiplicidad

Las ecuaciones más conocidas son las ecuaciones polinómicas. A ellas dedicamos esteapartado.

La resolución de las ecuaciones polinómicas es trivial cuando son de primer grado o desegundo grado. Pero la resolución de ecuaciones polinómicas de grado superior no es, por logeneral, inmediata. A continuación recordamos un teorema acerca de las raíces de polinomios�o soluciones de ecuaciones polinómicas, o ceros de funciones polinómicas: las tres cosasson lo mismo� con coe�cientes reales.

Teorema 2.1 Todo polinomio p(x) de grado n con coe�cientes reales �es decir, p(x) =

anxn + an−1x

n−1 + · · · + a1x + a0, con ai ∈ R para todo i = 0, 1, . . . , n y an = 0� puede

1En esta asignatura los términos polinomio y función polinómica serán prácticamente sinónimos.

26

expresarse de forma única (salvo el orden de los factores) de la siguiente forma, a la que se

llama descomposición canónica de p(x):

p(x) = an

u∏i=1

(x− ri)mi

v∏j=1

(x2 + sjx+ tj)nj . (2.2)

En esta descomposición cada factor x2 + sjx + tj es irreducible (sus raíces no son reales),

los coe�cientes ri, sj y tj son números reales para cada i = 1, 2, . . . , u y j = 1, 2, . . . , v, y los

exponentes mi y nj son números naturales tales que∑u

i=1 mi + 2∑v

j=1 nj = n.

Observe que el polinomio p(x) de la expresión (2.2) tiene como raíces reales a los númerosri (i = 1, 2, . . . , u), y como raíces imaginarias a las raíces de cada uno de los factores irre-ducibles x2 + sjx+ tj, a saber,

ρj =−sj − i

√4tj − s2j

2y ρj =

−sj + i√

4tj − s2j

2(j = 1, 2, . . . , v).

De ahí que se obtenga que

p(x) = an

u∏i=1

(x− ri)mi

v∏j=1

[(x− ρj)

nj(x− ρj)nj].

Al exponente mi de cada factor (x − ri) se le llama multiplicidad de la raíz ri. Del mismomodo, al exponente nj de los factores (x − ρj) y (x − ρj) se le llama multiplicidad de lasraíces ρj y ρj. En de�nitiva, se concluye el siguiente resultado.

Teorema 2.2 Si p(x) es el polinomio del teorema 2.1, entonces p(x) se puede descomponer

de forma única (salvo el orden de los factores) como producto de n polinomios de grado 1:

p(x) = an

w∏k=1

(x− zk)nk , donde zk ∈ C y nk ∈ N para todo k = 1, 2, . . . , w. (2.3)

Observe que en la expresión (2.3) se deber cumplir que∑w

k=1 nk = n.Con la ayuda del teorema 2.2, podemos de�nir de un modo más riguroso el concepto de

multiplicidad de una raíz de una ecuación polinómica, como se expresa a continuación.

De�nición 2.1 Un número α ∈ C es una raíz de multiplicidad k ∈ N de una ecuaciónpolinómica con coe�cientes reales p(x) = 0 si y solo si el polinomio p(x) puede descomponersede la siguiente forma:

p(x) = (x− α)kq(x), (2.4)

donde q(x) es un polinomio tal que q(α) = 0.

Como consecuencia del teorema 2.2 se puede a�rmar también que todo polinomio p(x)

de grado n tiene exactamente n raíces complejas (entre reales e imaginarias) si cada raíz secuenta tantas veces como indica su multiplicidad. Además, sus raíces imaginarias las tienea pares (si α es una raíz imaginaria de p(x), entonces α es también una raíz de p(x) con lamisma multiplicidad que α).

27

Si el polinomio es de coe�cientes enteros y tiene raíces enteras, el siguiente resultadoayuda a encontrar estas raíces y su multiplicidad. Con lo que se podrá descomponer elpolinomio como producto de factores utilizando la regla de Ru�ni, aunque lo más frecuenteserá que no se llegue hasta la descomposición canónica (2.2).

Teorema 2.3 Sea p(x) un polinomio de grado n ≥ 1 con coe�cientes enteros: p(x) = anxn+

an−1xn−1 + · · ·+ a1x+ a0, con ai ∈ Z (i = 0, 1, . . . , n) y an = 0. Entonces, las raíces enteras

de p(x), si existen, son divisores del término independiente a0.

Los siguientes ejemplos ilustran los tres resultados anteriores.

Ejemplo 2.1 Consideramos el polinomio

p(x) = 18 + 3x+ 23x2 − x3 − 27x4 − 14x5 − 2x6.

Aplicando el teorema 2.3, las raíces enteras de p(x), si existen, son divisores de 18. Losdivisores de 18 son los enteros ±1, ±2, ±3, ±6, ±9 y ±18. Probamos en primer lugar si losdivisores más pequeños son raíces de p(x). Así, aplicando Ru�ni se llega a que

p(x) = (x− 1)(x+ 2)(−9− 6x− 19x2 − 12x3 − 2x4).

El factor de grado 4 tiene a −9 como término independiente, cuyos divisores son ±1, ±3

y ±9. Pero a estas alturas ya se ha descartado que ±1 puedan ser raíces de dicho factor.Queda probar con ±3 y ±9, llegándose �nalmente a la siguiente descomposición

p(x) = −2(x− 1)(x+ 2)(x+ 3)2(x2 + 1/2),

que es la descomposición canónica que menciona el teorema 2.1. 2

Ejemplo 2.2 Aplicando el teorema 2.3 y la regla de Ru�ni al polinomio

p(x) = −2x5 − 8x4 + 2x3 + 12x2 − 36x,

es fácil obtener que su descomposición canónica es

p(x) = −2x(x+ 3)2(x2 − 2x+ 2).

El factor irreducible x2 − 2x+ 2 tiene como raíces imaginarias a 1 + i y 1− i. Luego

p(x) = −2x(3 + x)2(x− 1− i)(x− 1 + i).

es la descomposición de p(x) a la que se re�ere el teorema 2.2. 2

El teorema 2.3 permite hallar también las raíces enteras de los polinomios con coe�cientesracionales: basta multiplicar el polinomio por el mínimo común múltiplo de los denomi-nadores de sus coe�cientes, con lo que se obtiene un polinomio con coe�cientes enteros quetiene las mismas raíces que el polinomio de partida. Es lo que ocurre en el siguiente ejemplo.

28

Ejemplo 2.3 Dado el polinomio

p(x) =5

3− 65

18x− 305

9x2 − 105

6x3 − 5x4,

el mínimo común múltiplo de los denominadores de los coe�cientes es 18. Consideramos elpolinomio q(x) = 18p(x) = 30−65x−610x2−465x3−90x4 que, al ser de coe�cientes enteros,se le puede aplicar el teorema 2.3: las posibles raíces enteras de q(x) (y, por tanto, de p(x))serán divisores del número 30. Como en los ejemplos anteriores, se puede llegar sin di�cultada que q(x) = −90(x+ 2)(x+ 3)(x− 1

6)(x+ 1

3). Luego

p(x) = −5(x+ 2)(x+ 3)(x− 1

6

)(x+

1

3

)es la descomposición canónica de p(x). 2

Finalmente, presentamos otro resultado interesante acerca de las raíces de cualquier poli-nomio con coe�cientes reales.

Teorema 2.4 Sea p(x) = anxn + an−1x

n−1 + · · · + a1x + a0 un polinomio de grado n con

coe�cientes reales. Si p(x) tiene raíces reales, entonces el valor absoluto de cualquier raíz

real α de p(x) satisface que

|α| ≤ 1 +1

|an|max

0≤i≤n−1|ai|.

Por ejemplo, si se aplica el teorema 2.4 al polinomio del ejemplo 2.1, se obtiene que susraíces reales se encuentran en el intervalo [−1 − 27/2, 1 + 27/2]. De hecho vimos que esasraíces son 1, −2 y −3, que están en ese intervalo.

En la mayor parte de los casos, la obtención de las raíces de un polinomio con coe�cientesreales requiere la aplicación de métodos numéricos. Los que estudiaremos aquí serán útilespara muchos tipos de ecuaciones, no solo las polinómicas. En las secciones 2.3 y 2.4 desarro-llaremos algunos de estos métodos para aproximar las raíces reales de una ecuación (algunostambién son útiles para aproximar las raíces imaginarias, pero no nos detendremos en estepunto).

2.2.2. Raíces y su multiplicidad en ecuaciones no polinómicas

En este apartado estudiaremos cómo el concepto de multiplicidad de una raíz de unaecuación polinómica puede extenderse a ecuaciones no polinómicas.

En primer lugar, necesitamos introducir un resultado que caracteriza el concepto demultiplicidad de una raíz de una ecuación polinómica.

Teorema 2.5 Dados k ∈ N, α ∈ C y p(x) un polinomio de coe�cientes reales, entonces α

es una raíz de multiplicidad k de la ecuación p(x) = 0 si y solo si se satisface la siguiente

condición:

p(α) = p′(α) = · · · = p(k−1)(α) = 0 y p(k)(α) = 0.

La demostración de este teorema no es difícil pero, por brevedad, no la incluimos. Elteorema 2.5 puede utilizarse para extender la de�nición de multiplicidad de una raíz deuna ecuación polinómica a un tipo bastante más general de ecuaciones, como se establece acontinuación.

29

De�nición 2.2 Dados k ∈ N, α ∈ C y f una función k veces derivable en un entorno deα,2 se dice que α es una raíz de multiplicidad k de la ecuación f(x) = 0 si se cumple que

f(α) = f ′(α) = · · · = f (k−1)(α) = 0 y f (k)(α) = 0. (2.5)

Al igual que para polinomios, se puede hablar de raíces simples, dobles, triples. . . comotérminos sinónimos de los de raíces de multiplicidad 1, 2, 3. . . , respectivamente.

Ahora, la pregunta que lógicamente nos podemos hacer es si hay una caracterización dela de�nición de multiplicidad de una raíz que pueda asemejarse a la de�nición 2.1 (dada solopara ecuaciones polinómicas). La respuesta va a ser a�rmativa si exigimos alguna condiciónmás a la función f que de�ne la ecuación. Supongamos que α es una raíz de multiplicidad k

de la ecuación f(x) = 0, y que existe la derivada f (k+1) en un entorno Eα de α. Aplicando(2.5), el desarrollo de Taylor de orden k de la función f en el punto α (recuerde el teorema1.1) establece que

f(x) =f (k)(α)

k!(x− α)k +

f (k+1)(cx)

(k + 1)!(x− α)(k+1) = (x− α)k

(f (k)(α)

k!+

f (k+1)(cx)

(k + 1)!(x− α)

)para todo x ∈ Eα, donde cx es algún punto entre α y x. La igualdad obtenida es semejantea la igualdad (2.4): la función

q(x) =f (k)(α)

k!+

f (k+1)(cx)

(k + 1)!(x− α)

satisface que q(α) = f (k)(α)k!

= 0. La diferencia está en que la función q ya no es una funciónpolinómica (salvo que f lo fuese), puesto que cx = c(x) es una función desconocida de lavariable x.

A continuación mostramos un par de ejemplos de raíces de ecuaciones y sus multiplici-dades.

Ejemplo 2.4 Consideramos la función f(x) = x2 + ln(x + 1). Obsérvese que x = 0 es unaraíz de la ecuación f(x) = 0. Por otra parte, f ′(x) = 2x+1

x+1, con lo que f ′(0) = 1 = 0. Luego

x = 0 es una raíz simple de la ecuación f(x) = 0. 2

Ejemplo 2.5 Si g(x) = (x− 1)ex−1, es claro que x = 1 es una raíz de la ecuación g(x) = 0.Por otra parte, g′(x) = xex−1−1 y g′′(x) = (x+1)ex−1, con lo que g′(1) = 0 y g′′(1) = 1 = 0.Luego x = 1 es una raíz doble de la ecuación g(x) = 0. 2

En general, antes de resolver una ecuación se deberá determinar el número de solucionesque tiene �si es que tiene alguna� así como su localización (en un intervalo, si es real; en unaregión del plano, si es imaginaria); para todo lo cual será necesario utilizar algunos resultadosdel Análisis Real que no se incluyen en estos apuntes. En concreto, para obtener el número desoluciones reales de una ecuación del tipo (2.1) �sin determinar su multiplicidad� y sepa-rarlas en intervalos que contengan una sola de ellas, será necesario estudiar la continuidad, laderivabilidad y la monotonía de la función f , encontrar sus extremos locales, y quizá tambiénobtener algunos límites de esa función. Con la realización de algunos ejercicios de la sección2.6 se repasará este tipo de problemas.

2Por simplicidad, a lo largo de estos apuntes supondremos que un entorno Eα de un número real α es

cualquier intervalo de la forma (α− ε, α+ ε), con ε > 0.

30

2.3. Métodos de bisección, de la regula falsi y

de la secante

En esta sección desarrollamos tres métodos para calcular de modo aproximado las solu-ciones reales de la ecuación (2.1), bajo la hipótesis de que la función f que de�ne la ecuaciónes continua en el intervalo donde se ha localizado cada solución. Estos métodos no propor-cionan la multiplicidad de las soluciones aproximadas. Por otra parte, los dos primeros nopodrán aplicarse para aproximar raíces múltiples en las que la función f no cambie de signo.

Esos tres métodos, así como los que estudiaremos en las dos siguientes secciones, sonmétodos iterativos de resolución de ecuaciones (de sistemas de ecuaciones, en el caso de lasección 2.5), esto es, son métodos que generan una sucesión de números que se pretende queconverja a la solución buscada. Si esto ocurre, entonces se dice que el método iterativo esconvergente (convergente para el caso considerado: puede ocurrir que el mismo método parala misma solución sea convergente o no dependiendo del punto o puntos de partida donde seaplica el método).

Cuando nos re�ramos al error de la aproximación supondremos que se trata del errorabsoluto, que es el valor absoluto (o norma) de la diferencia entre la aproximación y lasolución de la ecuación (o del sistema de ecuaciones).

2.3.1. Método de bisección

El método iterativo más sencillo (pero también el más lento en cuanto a su convergenciaa la solución) es el conocido como método de bisección, que consiste en aplicar repetidamenteel teorema de Bolzano como se explica a continuación.

Dada la ecuación (2.1), supongamos que f es una función continua en el intervalo cerrado[a, b] y que f(a)f(b) < 0, de modo que puede asegurarse la existencia de al menos una soluciónde (2.1) en el interior del intervalo [a, b]. Para exponer este método y los demás que siguensupondremos, por simplicidad, que la ecuación (2.1) tiene una única solución s en el intervalo[a, b]. El método de bisección establece que una primera aproximación de s es el punto mediodel intervalo: x1 = (a + b)/2. Si evaluamos la función f en este punto nos encontramos contres posibilidades:

1. Que f(x1) = 0 (lo que ocurrirá raramente), en cuyo caso habremos hallado la solucións de manera exacta: s = x1.

2. Puede ocurrir que f(a)f(x1) < 0, en cuyo caso, aplicando de nuevo el teorema deBolzano, la solución s se encuentra en el interior del intervalo [a, x1].

3. La última posibilidad es que f(x1)f(b) < 0, lo que aseguraría que s se halla en elinterior del intervalo [x1, b].

En cualquiera de los dos últimos casos, llamamos [a2, b2] al intervalo que contiene lasolución y se repite el procedimiento: una segunda aproximación de s es entonces x2 =

(a2 + b2)/2, y volvemos a tener las tres mismas posibilidades que en el paso anterior. Laprimera es que s = x2; las otras dos llevarían a un nuevo intervalo [a3, b3] que contiene a s,y una nueva aproximación de s sería x3 = (a3 + b3)/2. Se razonaría de manera análoga paraobtener sucesivas aproximaciones x4, x5, etc. Luego el método de bisección �como todos

31

los métodos iterativos que estudiaremos en este capítulo y en el siguiente� producen unasucesión de números reales {xn}.3

En el método de bisección es claro que el error que se comete al aproximar la solución s

por xn, con n ∈ N, satisface que

|xn − s| < bn − an2

=bn−1 − an−1

22= · · · = b2 − a2

2n−1=

b− a

2n.

Con esta acotación se puede determinar fácilmente el número mínimo de iteraciones que esnecesario realizar para asegurar que el error de aproximación cometido sea menor que ciertacantidad positiva pre�jada ε, puesto que para que |xn − s| < ε es su�ciente exigir que

b− a

2n< ε lo que ocurre si y solo si n >

ln(b−aε

)ln 2

.

Por construcción �o lo que es lo mismo, como consecuencia de que |xn− s| ≤ (b− a)/2n

para todo n ∈ N�, es evidente que el método de bisección es siempre convergente, es decir,lımn→∞

xn = s.

A continuación se detalla un ejemplo de aplicación del método de bisección.

Ejemplo 2.6 Supongamos que se desea determinar una solución de la ecuación x3 = 3x−1

con un error de aproximación inferior a una milésima. Como la función f(x) = x3 − 3x + 1

es continua y se tiene que f(−2) = −1 < 0, f(0) = 1 > 0, f(1) = −1 < 0 y f(2) = 3 > 0,resulta que la ecuación dada, esto es, f(x) = 0, tiene tres soluciones reales situadas en losintervalos (−2, 0), (0, 1) y (1, 2).

A continuación aproximamos una de ellas por el método de bisección: por ejemplo, la quese encuentra en el intervalo (0, 1). Para asegurar que el error sea menor que una milésima setiene que cumplir que (1− 0)/2n < 10−3, es decir, n ≥ 10.

Luego los pasos de aplicación del método de bisección son en este caso los siguientes (enla última columna se expresa una cota del error de la aproximación obtenida en cada �la):

x1 = 1/2 f(x1) = −0,375 [a2, b2] = [0, 1/2] 1/2

x2 = 1/4 f(x2) = 0,265 . . . [a3, b3] = [1/4, 1/2] 1/4

x3 = 3/8 f(x3) = −0,072 . . . [a4, b4] = [1/4, 3/8] 1/8

x4 = 5/16 f(x4) = 0,093 . . . [a5, b5] = [5/16, 3/8] 1/16

x5 = 11/32 f(x5) = 0,009 . . . [a6, b6] = [11/32, 3/8] 1/32

x6 = 23/64 f(x6) = −0,031 . . . [a7, b7] = [11/32, 23/64] 1/64

x7 = 45/128 f(x7) = −0,011 . . . [a8, b8] = [11/32, 45/128] 1/128

x8 = 89/256 f(x8) = −0,0009 . . . [a9, b9] = [11/32, 89/256] 1/256

x9 = 177/512 f(x9) = 0,0042 . . . [a10, b10] = [177/512, 89/256] 1/512

x10 = 355/1024 1/1024

3En los métodos iterativos de resolución de sistemas obtendremos sucesiones de puntos de Rn. Por lo

general �tanto en el caso de sucesiones de números como en el de sucesiones de puntos de Rn�, conoceremos

el dominio de la sucesión. Si en alguna ocasión quisiéramos señalar expresamente cuál es ese dominio, lo

indicaríamos con una de las notaciones usuales de conjuntos: por ejemplo, la sucesión del método de bisección

se representaría por {xn : n ≥ 1}.

32

Con lo que se obtiene que x10 = 355/1024 = 0,3466796875 aproxima a una solución dex3 = 3x − 1 con un error inferior a una milésima (en realidad hemos visto que la cota esalgo inferior a una milésima, puesto que es 1/210 = 1/1024). De hecho, la solución exacta s

es (solo se muestran sus primeros decimales) 0,3472963553338606977 . . ., por lo que el errorcometido es aproximadamente igual a 6 diezmilésimas: |x10 − s| ≃ 0,0006166678. 2

A la hora de programar el método de bisección se podrá exigir como criterio de parada

(esto es, como criterio que determina el número de iteraciones a realizar) el ya comentado(que el error de aproximación sea menor que cierta cantidad) o también que la cantidad|f(xn)| sea menor que cierta cantidad positiva δ. Para detener la iteración se puede exigiruna u otra condición de parada, o incluso puede ser conveniente exigir que se cumplan ambas.El segundo criterio de parada mencionado podrá utilizarse también en todos los métodos queestudiaremos en este capítulo.

Como se ha mostrado, el método de bisección es convergente, pero esto no signi�ca quetoda aproximación dada por el método sea mejor que las anteriores. De hecho, en el ejemplo2.6 ocurre que |x8 − s| < |x10 − s|, como puede comprobarse utilizando Mathematica (perosi se calcula x11 en dicho ejemplo, puede comprobarse igualmente que x11 aproxima a lasolución de la ecuación mejor que todas las aproximaciones anteriores).

Conviene resaltar que el teorema de Bolzano puede utilizarse para aproximar raíces deecuaciones sin necesidad de aplicar el método de bisección: pueden elegirse las sucesivasaproximaciones sin que tengan que ser obligatoriamente el punto medio del intervalo dondese ha localizado la solución. Esto puede ser útil para ejercicios que necesiten pocas iteracionesy estas se calculen sin la ayuda del ordenador, como ocurre en el siguiente ejemplo.

Ejemplo 2.7 Supongamos que se desea determinar una solución negativa de la ecuación

5 +3

x

[1−

(100

100− x

)360]= 0,

con dos dígitos decimales correctos. Como la función f(x) = 5+ 3x

[1−

(100

100−x

)360]es continua

en R−, y se tiene que f(−1) ≃ 2,08 > 0 y f(−0,5) ≃ −0,0037 < 0, resulta que la ecuacióntiene al menos una solución en el intervalo (−1,−0,5). Como f(−0,5) es un valor bastantepróximo a 0 se puede sospechar que x = −0,5 puede estar próximo a una solución de laecuación (no es seguro que sea así, con lo datos que tenemos hasta ahora, pero �por si lofuese� vale la pena intentar buscar una solución cerca de x = −0,5).

Aquí el método de bisección impondría tomar x1 = −0,75, pero los cálculos del párrafoanterior sugieren que la primera aproximación x1 a elegir en el intervalo (−1,−0,5) se sitúemucho más próxima al extremo derecho del intervalo, como por ejemplo x1 = −0,55. Comof(x1) = f(−0,55) ≃ 0,3 > 0, hay una solución s el intervalo (−0,55,−0,5), y parece probableque esté más próxima a x = −0,5 que a x1. Tomamos entonces x2 = −0,51. Como f(x2) =

f(−0,51) ≃ 0,06 > 0, hay una solución s el intervalo (−0,51,−0,50), y por tanto s es dela forma s = −0,50 . . ., es decir, una aproximación de s con dos decimales correctos ess∗ = −0,50. Observe que si se hubiera usado el método de bisección habríamos necesitadoobtener un número mayor aproximaciones intermedias. 2

33

2.3.2. Método de la regula falsi

El método de la regula falsi requiere las mismas hipótesis que el método de bisección ytambién aplica repetidamente el teorema de Bolzano. Será igualmente convergente, pero suconvergencia será en general más rápida que la del método de bisección.

Supongamos de nuevo que la función f de la ecuación (2.1) es continua en un intervalocerrado [a, b], que f(a)f(b) < 0 y que f tiene un único cero s en (a, b). El método de laregula falsi establece que una primera aproximación x1 de s es el punto de corte de la rectaque une los dos puntos extremos de la grá�ca de f en [a, b] con el eje de abscisas. Como laecuación de dicha recta es

y − f(a)

f(b)− f(a)=

x− a

b− a,

es inmediato obtener que4

x1 = a− f(a)(b− a)

f(b)− f(a).

Si evaluamos la función f en este punto hay tres posibilidades idénticas a las que se con-sideraron en el método de bisección. Entonces, si x1 = s, se denota [a2, b2] al intervalo quecontiene la solución ([a, x1] o [x1, b]) y se repite el procedimiento: una segunda aproximaciónde s es entonces

x2 = a2 −f(a2)(b2 − a2)

f(b2)− f(a2),

y volvemos a tener las tres mismas posibilidades que en el paso anterior. De manera análogase obtendrían las sucesivas aproximaciones. Por tanto, la iteración del método de la regula

falsi puede expresarse como sigue:

xn = an −f(an)(bn − an)

f(bn)− f(an), n ∈ N, (2.6)

donde a1 = a y b1 = b, mientras que ai y bi (para todo i > 1) se obtienen como hemosexplicado.

El estudio del error de aproximación del método de la regula falsi y de los demás méto-dos que estudiaremos más adelante en este capítulo es complicado, por lo que lo omitiremos.Pero al menos se puede a�rmar que el error de aproximación del método de la regula falsi

tiende a cero a medida que aumentan las iteraciones, es decir, el método de la regula falsi

es convergente. Además, por lo general converge con bastante más rapidez que el método debisección. El criterio de parada suele ser que el valor absoluto de la diferencia de dos aproxi-maciones consecutivas (cantidad que suele ser menor que el error de la última aproximaciónobtenida) sea menor que cierta cantidad ε. O también podría ser el ya mencionado de quela cantidad |f(xn)| sea menor que cierta cantidad positiva δ. También puede exigirse quese cumplan ambos criterios antes de parar el algoritmo. Ambos criterios de parada seránválidos también para los demás métodos de aproximación de soluciones que estudiaremos eneste capítulo.

4Otras fórmulas equivalentes son x1 = b − f(b)(b− a)

f(b)− f(a)y x1 =

af(b)− bf(a)

f(b)− f(a), pero esta última suele

producir más errores de redondeo que las otras dos.

34

En general (en este y en los demás métodos), si se utiliza un criterio de parada puede serrecomendable comprobar que la condición del otro criterio se cumple satisfactoriamente. Porejemplo, si se utiliza el primer criterio de parada expuesto en el método de la regula falsi, sepuede obtener el valor |f(xn)| (para la última aproximación calculada xn) y comprobar quees su�cientemente pequeño conforme a la precisión que se necesita.

A continuación ilustramos el método de la regula falsi aplicándolo a la misma ecuacióndel ejemplo 2.6.

Ejemplo 2.8 Dada la ecuación x3 = 3x − 1, equivalentemente f(x) = x3 − 3x + 1 = 0,aplicamos el método de la regula falsi para aproximar la solución de esa ecuación que seencuentra en el intervalo (0, 1). Establecemos como criterio de parada que la diferencia dedos aproximaciones consecutivas sea menor que 10−6. Los pasos de aplicación de este métodoproducen los siguientes resultados: el intervalo de partida es [a1, b1] = [0, 1]. Por tanto, paran = 1 la igualdad (2.6) se traduce en que

x1 = 0− f(0)(1− 0)

f(1)− f(0)=

1

2.

Como f(x1) = −0,375 < 0, la solución se encuentra en el intervalo [a2, b2] = [a1, x1] = [0, 1/2].Esta y las siguientes iteraciones se describen resumidamente a continuación:

x1 = 1/2 = 0,5 f(x1) = −0,375 [a2, b2] = [a1, x1]

x2 = 4/11 ≃ 0,3636 |x2 − x1| ≃ 1,36 · 10−1 f(x2) ≃ −4,28 · 10−2 [a3, b3] = [a1, x2]

x3 = 121/347 ≃ 0,3487 |x3 − x2| ≃ 1,49 · 10−2 f(x3) ≃ −3,71 · 10−3 [a4, b4] = [a1, x3]

x4 ≃ 0,347414 |x4 − x3| ≃ 1,29 · 10−3 f(x4) ≃ −3,12 · 10−4 [a5, b5] = [a1, x4]

x5 ≃ 0,347306 |x5 − x4| ≃ 1,08 · 10−4 f(x5) ≃ −2,61 · 10−5 [a6, b6] = [a1, x5]

x6 ≃ 0,34729718 |x6 − x5| ≃ 9,07 · 10−6 f(x6) ≃ −2,19 · 10−6 [a7, b7] = [a1, x6]

x7 ≃ 0,34729642 |x7 − x6| ≃ 7,60 · 10−7 f(x7) ≃ −1,83 · 10−7 [a8, b8] = [a1, x7]

Con la ayuda de Mathematica puede probarse que el error que se produce al tomar x7

como aproximación de la solución es del orden de 6,95 · 10−8: observe que es unas diez vecesmenor que |x7 − x6|. 2

El hecho de que el valor de ai se haya mantenido constante desde el principio en el ejemplo2.8, no quiere decir que esto suceda siempre en el método de la regula falsi. Aunque sí sueleocurrir que, a partir de cierto paso, los valores de uno de los extremos de los intervalos [ai, bi]permanezca invariable.

2.3.3. Método de la secante

El método de la secante se basa en una idea similar a la del método de la regula falsi

pero, a diferencia de este y del método de bisección, no se parte necesariamente en cadapaso de un intervalo que contenga a la solución, es decir, no utiliza el teorema de Bolzano.Como inconveniente tiene que no siempre es convergente; sin embargo, cuando lo es, suconvergencia es en general bastante más rápida que la del método de la regula falsi una vezrealizadas unas pocas iteraciones.

35

A continuación se detalla en qué consiste el método de la secante. Supongamos que lafunción f de la ecuación (2.1) es continua en un intervalo I donde tiene una solución s quese desea aproximar. Supongamos que x−1 y x0 son dos aproximaciones de s (suele ocurrirque f(x−1)f(x0) < 0 y que la solución s que se busca está entre esas dos aproximaciones,pero no siempre será así) tales que f(x−1) = f(x0). En este caso, la recta que contiene a lospuntos (x−1, f(x−1)) y (x0, f(x0)) corta al eje de abscisas en un punto x1, que se pretendeque sea una nueva aproximación de s.

El método de la secante repite el proceso, si es posible, con los puntos (x0, f(x0)) y(x1, f(x1)), de manera que x2 sería el punto de corte de la recta que pasa por los puntos(x0, f(x0)) y (x1, f(x1)) con el eje de abscisas; y así sucesivamente. Por tanto, las sucesivasaproximaciones se obtienen como en el método de la regula falsi, aunque ahora la notaciónse simpli�ca:

xn+1 = xn −f(xn)(xn − xn−1)

f(xn)− f(xn−1), n = 0, 1, 2, . . . 5

Los criterios de parada del método de la secante son los mismos que los del método de laregula falsi. La diferencia está en que cuando el método no converja el algoritmo de la secantecaerá en un bucle in�nito o, más frecuentemente, producirá un error: porque la sucesión delos puntos {xn} diverja, o porque algún xn caiga fuera del intervalo I donde f es continua,o porque antes de converger suceda que f(xn) = f(xn−1) para algún valor de n.

En el siguiente ejemplo se aproxima la solución de la ecuación de los ejemplos 2.6 y 2.8utilizando el método de la secante. En este caso la solución s que se desea aproximar ha sidolocalizada en un intervalo [a, b] tal que f(a)f(b) < 0, y puede considerarse que los extremosdel intervalo son dos aproximaciones de partida de la solución: x−1 = a y x0 = b (en caso deque el método no converja con esta elección, se puede probar si converge tomando x−1 = b

y x0 = a). Así, en este caso particular que estamos considerando, sucede que las primeraaproximación x1 del método de la secante coincide con primera aproximación del método dela regula falsi ; teóricamente podrían coincidir algunas aproximaciones más (x1, x2,. . . ) peronormalmente serán muy pocas.

Ejemplo 2.9 Aplicamos el método de la secante para aproximar la solución de la ecuaciónx3 = 3x−1 que se encuentra en el intervalo [0, 1]. Al igual que en el ejemplo 2.8, se establececomo criterio de parada que la diferencia de dos aproximaciones consecutivas sea menor que10−6. Partiendo de x−1 = 0 y x0 = 1, las primeras iteraciones que produce el método de lasecante aplicado a esa ecuación son las siguientes:

x1 = 1/2 = 0,5 |x1 − x0| = 5 · 10−1 f(x1) = −3,75 · 10−1

x2 = 1/5 = 0,2 |x2 − x1| = 3 · 10−1 f(x2) = 4,08 · 10−1

x3 = 31/87 ≃ 0,3563 |x3 − x2| ≃ 1,56 · 10−1 f(x3) ≃ −2,37 · 10−2

x4 ≃ 0,347731 |x4 − x3| ≃ 8,59 · 10−3 f(x4) ≃ −1,15 · 10−3

x5 ≃ 0,34729478 |x5 − x4| ≃ 4,37 · 10−4 f(x5) ≃ 4,14 · 10−6

x6 ≃ 0,3472963556 |x6 − x5| ≃ 1,57 · 10−6 f(x6) ≃ −7,12 · 10−10

x7 ≃ 0,3472963553338609 |x7 − x6| ≃ 2,70 · 10−10 f(x7) ≃ −4,44 · 10−16

5Una fórmula equivalente es: xn+1 = xn−1 −f(xn−1)(xn − xn−1)

f(xn)− f(xn−1), n = 0, 1, 2, . . .

36

Con la ayuda de Mathematica puede probarse que el error que se produce al tomar x7

como aproximación de la solución es del orden de 2,023 · 10−16: observe que es un millón deveces menor que |x7−x6| y unas 108 más pequeño que el obtenido con la séptima aproximacióndel método de la regula falsi (obtenida en el ejemplo 2.8). Es más, en los primeros pasos deeste ejemplo el método de la secante obtuvo peores resultados que el método de la regula falsi,pero la �aceleración� de la convergencia del método de la secante ha sido bastante mayor.2

2.4. Método de Newton-Raphson

En esta sección desarrollamos un nuevo método iterativo para calcular de modo aproxima-do las soluciones de la ecuación (2.1): el método de Newton-Raphson. Una primera ventajade este método respecto de los estudiados en la sección 2.3 es que cada aproximación seobtiene solo a partir de la anterior (en los otros métodos se usaban dos aproximaciones an-teriores). Pero lo más importante es que, cuando el método sea convergente, su velocidad deconvergencia será casi siempre bastante mayor que en dichos métodos. A diferencia de estos,la aplicación del método de Newton-Raphson requerirá una hipótesis más sobre la funciónf , que deberá ser derivable; y las iteraciones del método requerirán evaluar no solo f sinotambién su derivada.

El método de Newton-Raphson tiene como los anteriores una clara interpretación geo-métrica. Se puede describir como sigue. Si se desea aproximar una determinada solución s

de la ecuación (2.1) por el método de Newton-Raphson, se deben seguir los siguientes pasos.

1. Se parte de un punto x0 que esté más o menos próximo a la solución s (la necesidadde una mayor o menor proximidad a s depende de la función considerada y �cuandola ecuación tiene más de una solución� de la solución a aproximar). Frecuentementese parte de una situación en la que se ha localizado s en un intervalo que contienea una única solución de la ecuación, en cuyo caso el punto x0 se puede elegir entrelos del intervalo (a veces se elige el punto medio del intervalo, buscando un error departida menor que la mitad de la longitud del intervalo; pero hay otras motivacionesque pueden in�uir en la elección de este punto inicial).

2. Se traza la recta tangente a la grá�ca de f en el punto (x0, f(x0)) (supuesto que existay no sea horizontal). Se pretende que la intersección de esta recta con el eje de abscisasproporcione una primera aproximación x1 de s.

3. Sobre x1 se repite el proceso realizado con x0, obteniéndose x2; de igual modo seobtienen sucesivamente x3, x4, etc. Es fácil probar que el método descrito produce unasucesión de números que se obtiene mediante la siguiente fórmula:

xn+1 = xn −f(xn)

f ′(xn), n = 0, 1, 2, . . . (2.7)

4. Si se observa que la sucesión {xn} converge a s, se realiza la iteración (2.7) hastaque se obtenga una aproximación con la exactitud deseada (a medida que avanza laiteración, los dos últimos valores de la sucesión van coincidiendo en un mayor númerode dígitos, que pueden considerarse correctos). Pero si resulta que {xn} no converge

37

a s pero converge a s ∈ R, entonces s es otra solución de la ecuación (2.1). Tantosi sucede esto último como si {xn} no converge a ningún número real, habría quebuscar un punto de partida x0 distinto del anterior (en principio, más próximo a s;por ejemplo, el nuevo valor de x0 puede obtenerse a partir del su valor anterior y deotro punto tras realizar una o más iteraciones de alguno de los tres métodos estudiadosen la sección anterior) de modo que la sucesión (2.7) sí converja a s. Casi siempreserá posible encontrar ese valor adecuado de partida x0; pero solo estudiaremos cómohacerlo en ejercicios concretos �a mano y con Mathematica�, sin introducir nuevosteoremas que alargarían demasiado esta exposición.

Los criterios de parada que utilizaremos para lograr un error adecuado a nuestros deseosserán los mismos que los considerados para los métodos de la secante y de la regula falsi.

A continuación aplicamos el método de Newton-Raphson a la ecuación estudiada en losejemplos 2.6, 2.8 y 2.9.

Ejemplo 2.10 Se desea aplicar el método de Newton-Raphson para aproximar la soluciónde la ecuación x3 = 3x − 1, equivalentemente f(x) = x3 − 3x + 1 = 0, que se encuentra enel intervalo (0, 1). Como en los ejemplos 2.8 y 2.9, se establece como criterio de parada quela diferencia de dos aproximaciones consecutivas sea menor que 10−6. Tomando x0 = 1/2 seobtienen los siguientes resultados:

x0 = 1/2 = 0,5 f(x0) = −3,75 · 10−1

x1 = 1/3 ≃ 0,333 |x1 − x0| = 1,66 · 10−1 f(x1) = 3,70 · 10−2

x2 = 25/72 ≃ 0,347222 |x2 − x1| ≃ 1,39 · 10−2 f(x2) ≃ 1,96 · 10−4

x3 ≃ 0,34729635316 |x3 − x2| ≃ 7,41 · 10−5 f(x3) ≃ 5,72 · 10−9

x4 ≃ 0,3472963553338607 |x4 − x3| ≃ 2,17 · 10−9 f(x4) ≃ −2,22 · 10−16

Con la ayuda de Mathematica puede probarse que el error absoluto que se produce altomar x4 como aproximación de la solución es aproximadamente igual a 2,3 · 10−18: observeque es casi mil millones de veces menor que |x4 − x3| y es menor también que el obtenidocon la séptima aproximación del método de la secante. 2

El siguiente ejemplo obtiene la fórmula de Herón (año 100 a. C.) para aproximar elcálculo de la raíz cuadrada de un número mediante sumas, restas, productos y divisiones.En su tiempo no se conocía el algoritmo para la obtención de la raíz cuadrada.

Ejemplo 2.11 Se aplica el método de Newton-Raphson para aproximar la raíz cuadradapositiva de un número positivo a, esto es, la solución positiva de la ecuación x2−a = 0. Herónpartía del mayor número natural x0 tal que x2

0 < a. Al aplicar el método de Newton-Raphsona esa ecuación se obtiene la fórmula de Herón:

xn+1 = xn −x2n − a

2xn

=x2n + a

2xn

, n = 0, 1, 2, . . .

Por ejemplo, para a = 7 esta fórmula produce las siguientes aproximaciones de√7:

x0 = 2, x1 =11

4= 2,75, x2 =

233

88≃ 2,6477, x3 =

108497

41008≃ 2,645752.

La última aproximación satisface que x3 −√7 ≃ 7,37 · 10−7. 2

38

Método de Newton-Raphson modi�cado

Si la multiplicidad r de la raíz que se desea obtener aproximar es mayor que 1, entonceshay una variación del método de Newton-Raphson que converge con más rapidez a esa raízmúltiple. Consiste en multiplicar el sustraendo de la fórmula (2.7) por la multiplicidad dedicha raíz, es decir,

xn+1 = xn − rf(xn)

f ′(xn), n = 0, 1, 2, . . . (2.8)

Es frecuente que, a priori, no se conozca dicha multiplicidad r. Si al aplicar el método(2.7) observamos que converge lentamente a la raíz buscada s, entonces es conveniente probarel método (2.8) para r = 2. Si este método convergiera aún más lentamente a s (o no con-vergiera) es muy probable que la raíz s sea simple. Pero lo más frecuente será que el método(2.8) converja con más rapidez, en cuyo caso casi podríamos asegurar que la multiplicidadde la raíz s es al menos 2. Entonces deberíamos utilizar de nuevo el método (2.8) pero conr = 3, que puede converger a s más lentamente que con r = 2 (o incluso no converger), loque haría sospechar que la multiplicidad de la raíz s es 2, o puede converger a s con másrapidez, en cuyo caso la multiplicidad de la raíz s es probablemente 3 o mayor que 3. Elproceso continúa hasta que encontremos un valor de r que converja más lentamente a s quecon el valor anterior de r (o incluso que no converja).

Para estar más seguros de las conclusiones obtenidas en el párrafo anterior, si se partede distintos puntos x0 y para la gran mayoría de ellos el método (2.8) converge más rápidopara cierto valor r = r0, es bastante seguro que la multiplicidad de la raíz s sea r0.

2.5. Método de Newton-Raphson para sistemas

de ecuaciones no lineales

En esta sección introducimos un método iterativo para calcular de modo aproximadolas soluciones de sistemas de ecuaciones no lineales con igual número de ecuaciones quede incógnitas.6 El problema que nos encontramos en este caso es que la determinación delnúmero de soluciones así como su localización es en general mucho más complicada quecuando se tiene una sola ecuación con una incógnita: de hecho, no abordaremos aquí eseproblema.

Para facilitar la comprensión y la expresión del problema, consideraremos primero el casomás simple (dos ecuaciones con dos incógnitas), en cuyo caso el sistema puede expresarse dela siguiente forma:

f(x, y) = 0,

g(x, y) = 0.

}(2.9)

6Los sistemas de ecuaciones lineales se saben resolver de manera exacta, por ejemplo por el método de

Gauss. De todos modos, algunos problemas numéricos relacionados con este tipo de sistemas se estudiarán

en el capítulo 3.

39

Para lo que sigue, se necesitará la siguiente función matricial:

J(x, y) =

∂f

∂x(x, y)

∂f

∂y(x, y)

∂g

∂x(x, y)

∂g

∂y(x, y)

. (2.10)

A J(x, y) se le llama matriz jacobiana de la función vectorial (f, g) en el punto (x, y).El método de Newton-Raphson se describe en el siguiente teorema, que asegura su conver-

gencia. En su enunciado se exige que las funciones f y g sean de clase C2 en cierto dominio,lo que signi�ca que todas las derivadas parciales de primer y segundo orden de esas funcionesexisten y son continuas en dicho dominio.

Teorema 2.6 Sean f y g dos funciones de clase C2 en cierto dominio D que contiene en su

interior una solución (s, t) del sistema (2.9). Si det(J(s, t)) = 0, entonces existe un entorno

E del punto (s, t) tal que para cualquier (x0, y0) ∈ E la sucesión de puntos obtenida mediante

la siguiente iteración(xn+1

yn+1

)=

(xn

yn

)− J−1(xn, yn)

(f(xn, yn)

g(xn, yn)

), n = 0, 1, 2, . . . (2.11)

converge a (s, t).

Como (s, t) es desconocido, para comprobar la condición det(J(s, t)) = 0 será necesarioprobar que det(J(x, y)) = 0 en cierto dominio de R2 que contenga en su interior a (s, t). Elteorema tampoco proporciona ninguna información acerca del entorno E. En la práctica, unavez seleccionado el punto (x0, y0), si la sucesión {(xn, yn)} obtenida mediante (2.11) convergea un punto del dominio D, este punto es la solución buscada; si no ocurriera esto, habríaque encontrar otro punto de partida (x0, y0) que se aproxime su�cientemente a la solución.

En el caso general de un sistema m×m, con m ≥ 2, a saber,

F1(x1, x2, . . . , xm) = 0,

F2(x1, x2, . . . , xm) = 0,...

...

Fm(x1, x2, . . . , xm) = 0,

(2.12)

se veri�ca un resultado análogo al teorema 2.6: si las funciones Fi (i = 1, 2, . . . ,m) son declase C2 en cierto dominio D que contiene en su interior una solución s = (s1, s2, . . . , sm)

del sistema (2.12), y la función matricial

J(x1, . . . , xm) =

∂F1

∂x1

(x1, . . . , xm)∂F1

∂x2

(x1, . . . , xm) · · · ∂F1

∂xm

(x1, . . . , xm)

∂F2

∂x1

(x1, . . . , xm)∂F2

∂x2

(x1, . . . , xm) · · · ∂F2

∂xm

(x1, . . . , xm)

......

...

∂Fm

∂x1

(x1, . . . , xm)∂Fm

∂x2

(x1, . . . , xm) · · · ∂Fm

∂xm

(x1, . . . , xm)

40

satisface que det(J(s)) = 0, entonces existe un entorno E del punto s tal que, para cualquierx(0) = (x

(0)1 , x

(0)2 , . . . x

(0)m ) ∈ E, la sucesión de puntos obtenida mediante la iteración

x(n+1) = x(n) − J−1(x(n)) · F (x(n)), n = 0, 1, 2, . . . , (2.13)

donde F = (F1, F2, . . . Fm), converge a s.Los comentarios hechos tras el teorema 2.6 son igualmente válidos para el caso m > 2.Cuandom sea bastante mayor que 2, normalmente el cálculo de la inversa será complicado

y puede producir errores de redondeo. En este caso la igualdad (2.13) no será la mejor formade obtener la sucesión de aproximaciones x(n). Para evitar el cálculo de la inversa puedemultiplicarse por J(x(n)) en la igualdad (2.13), obteniéndose

J(x(n)) · (x(n+1) − x(n)) = −F (x(n)), n = 0, 1, 2, . . . (2.14)

Llamando u(n) = x(n+1) − x(n), la ecuación (2.14) se transforma en el sistema lineal m×m

J(x(n)) · u(n) = −F (x(n)), n = 0, 1, 2, . . . , (2.15)

cuyas incógnitas son lasm componentes del vector u(n). Una vez resuelto el sistema (2.15) poralgún método exacto (como el de Gauss) o incluso mediante un método numérico (como se hadicho, los métodos numéricos de resolución de sistemas lineales se estudiarán en el capítulo3), la aproximación x(n+1) de (2.12) se obtiene sumándole a la anterior aproximación x(n) elvector u(n) recién calculado:

x(n+1) = x(n) + u(n), n = 0, 1, 2, . . .

Los criterios de parada para lograr un error adecuado a nuestros deseos serán semejantesa los considerados en la sección anterior, sustituyendo los valores absolutos por normas.En concreto, el primer criterio de parada consiste en que la norma de la diferencia de dosaproximaciones consecutivas sea su�cientemente pequeña (∥x(n+1) − x(n)∥ < ε, para ciertovalor pre�jado ε). Y el segundo criterio consiste en que la norma de la función que de�neel sistema, aplicada en una aproximación x(n), sea su�cientemente pequeña (∥F (x(n))∥ < δ,para cierto δ).

2.6. Ejercicios

1. Determine para qué números reales no nulos c y d la ecuación x2 − 4πx = c + d cos x

tiene una raíz de multiplicidad 4. ¾Cuál es esa raíz?

2. a) ¾Para qué valores k ∈ R la ecuación 2 cos 2x+ 4x+ k = 0 tiene una raíz triple?

b) ¾Para qué valores k ∈ R la raíz x = 0 de la ecuación ex−1 sen 2x = kx es doble?

3. Una persona necesita un préstamo de 150000 euros para la compra de un piso, préstamoque desea amortizar mediante pagos mensuales durante 30 años. Sea i el interés mensualdel préstamo (en tanto por ciento, expresado con no más de dos dígitos decimales). Lafórmula de amortización del préstamo es la siguiente:

150000 =100M

i

[1−

(100

100 + i

)360],

41

donde M es la cantidad mensual a pagar. Si dicha persona no está dispuesta a pagarmás de 900 euros mensuales, ¾cuál es el máximo interés mensual que puede asumir?

4. Pruebe que la ecuación 6x3 − 5x2 + 7x = 2 tiene una única solución. Aproxime esasolución con un error menor que 10−2 utilizando el método de bisección.

5. Calcule de manera exacta las raíces del polinomio p(x) = −0,4x2+2,2x+4,7 utilizandola fórmula usual. Realice tres iteraciones del método de bisección para determinar laraíz mayor a partir de un intervalo de amplitud 10−1 que contenga dicha raíz. Obtengauna cota del error y compárela con el verdadero error después de cada iteración. Observesi ese error disminuye en cada aplicación del método o no.

6. Haga una representación aproximada de la función f(x) = x sen x que contenga variosde sus ceros. ¾Qué tiene de especial el punto x = 0 si lo comparamos con los demásceros de f? ¾Podría utilizarse el método de bisección para aproximar x = 0 comosolución de la ecuación x senx = 0?

7. Localice todas las soluciones de la ecuación 1/x = tg x en intervalos que contengan unasola. ¾Se podría aplicar el método de bisección para aproximarlas? Aproxime algunade las soluciones mediante el método de bisección de manera que el error cometido seamenor que 10−2. Repita el ejercicio para la ecuación 1/x = 2x.

8. Utilice el método de bisección para encontrar una aproximación, con un error máximode una centésima, de la abscisa del punto donde se cortan las grá�cas de las funcionesf(x) =

√x2 + 1 y g(x) = tg x, de�nidas en el intervalo (0, π/2).

9. [Mathematica] Utilizando el método de bisección, aproxime con un error menor que10−5 una solución de la ecuación

√x− 3 = 3 ln(x − 10). Programe el algoritmo de

modo que: (a) se visualicen las sucesivas aproximaciones; (b) la aproximación �nal serepresente con al menos 7 dígitos de mantisa; (c) junto a cada aproximación se muestreuna cota del error, de manera que esta cota sea menor que 10−5 solo en el último paso.

10. Se considera la ecuación lnx = ax+ b, donde a y b son dos números reales.

a) Determine si es posible que esta ecuación tenga una raíz triple en algún caso, paraciertos valores de a y b.

b) Si b = −2, encuentre para qué valor o valores de a ocurre que la ecuación tieneuna raíz doble. En este caso, determine el valor exacto de dicha raíz.

c) Si a = 1 y b = −2, encuentre un intervalo de la forma [n, n + 1], con n ≥ 1,que contenga una raíz de la ecuación. Aplique el método de la �regula falsi� paraaproximar esta raíz hasta que la diferencia de dos aproximaciones consecutivassea menor que 10−3.

11. Aproxime las raíces de las siguientes ecuaciones utilizando los métodos de la regula-falsiy de la secante, con al menos tres iteraciones cada uno, aplicándolos al intervalo quele parezca conveniente:

(a) ex = 2− senx (b) ex = tg x (c) x3 − 12x2 + 16x+ 1 = 0.

42

12. Una bala se disparó verticalmente en el aire y está descendiendo a su velocidad terminalv, que es la dada por la expresión 1,4 · 10−2v1,5 + 21,15v = v2. Determine aproximada-mente v por el método de la secante si una estimación burda de v es v = 20 m/s.

13. a) Si f(x) = x3 − x − 3, pruebe que la ecuación f(x) = 0 tiene una única soluciónreal.

b) [Mathematica] Aproxime esta solución mediante el método de la secante, concondición de parada |f(x)| < 10−6. ¾Cuál es la solución aproximada que se obtienesi el intervalo de partida es [−10, 20]? ¾Y si es [1, 2]? ¾Qué resultado proporcionanlas tres instrucciones directas de Mathematica? Comente los resultados.

14. a) Demuestre que la ecuación x− cos x = 3 tiene una única solución real.

b) [Mathematica] Aproxime esta solución mediante el método de la secante, concondición de parada |f(x)| < 10−6.

15. [Mathematica] Considere la ecuación x3 = x cos x+20√x4 + x2 e implemente el método

de la secante con la condición de parada |xn+1 − xn| < 10−6 para aproximar la únicasolución no nula de dicha ecuación.

16. [Mathematica] La ecuación x3 + ex = 10 tiene una única raíz real. Determine unintervalo en el cuál esté la raíz y utilice los métodos de bisección, de la regula falsi,de la secante y de Newton-Raphson hasta que la diferencia entre dos aproximacionesconsecutivas sea menor que 10−4. Comente los resultados.

16b. [Mathematica] Utilice los métodos de bisección, regula falsi, secante y Newton-Raphsonpara obtener aproximadamente las soluciones de la ecuación x3 − 10 = ex, de maneraque la diferencia entre dos aproximaciones consecutivas sea menor que 10−4. Comentelos resultados.

17. Mediante el método de Newton-Raphson, aproxime las raíces más cercanas a 4,5 y a7,7 de la ecuación x = tg x realizando 3 y 5 iteraciones respectivamente. Sugiera unaposible cota del error cometido para cada aproximación. Obtenga en ambos casos elvalor absoluto de la diferencia de las dos últimas aproximaciones.

18. Pruebe que la ecuación x2 + ln x = 0 tiene una única solución. Aproxime la soluciónrealizando tres iteraciones del método de Newton-Raphson. Dé una opinión razonadaacerca de la exactitud de la aproximación alcanzada.

19. Se considera la ecuación x = cosx + 3. En el ejercicio 14 se demostró que tiene unaúnica solución real. (a) Localice esa solución en un intervalo de la forma [n, n+1], donden es un número entero. (b) Utilice el método de Newton-Raphson para aproximarlacon al menos dos decimales de exactitud. Indicación: habrá que obtener xk y xk+1 demodo que |xk+1 − xk| = 0,00 . . .

20. [Mathematica] Obtenga mediante el método de Newton-Raphson una aproximación decada una de las raíces de la ecuación xe−x−x2+1 = 0. Cada aproximación debe tenerun mínimo de 5 dígitos correctos.

43

21. Encuentre una ecuación f(x) = 0 tal que una de sus soluciones sea√a (a>0) y de modo

que, aplicando el método de Newton-Raphson a esa ecuación,√a pueda aproximarse

mediante operaciones aritméticas básicas (sumas, restas, multiplicaciones y divisiones).En particular, obtenga un valor aproximado z de

√26 tal que |s2 − 26| < 10−4.

22. Planteando la ecuación adecuada, el método de Newton-Raphson permite aproximarel inverso de un número real no nulo sin efectuar divisiones. En particular, calcule unaaproximación de π−1 mediante 3 iteraciones partiendo de x0 = 0,3 y compárela con elvalor de π−1 que le proporciona su calculadora u ordenador. Indicación para el casode que no encuentre esa ecuación: tome 1

ax= 1.

23. La velocidad de caída v (en m/s) de un paracaidista se obtiene mediante la siguienteigualdad:

v =gm

c

(1− e−ct/m

),

donde t es el tiempo transcurrido desde que se lanzó al vacío (en segundos), m es lamasa del paracaidista (en kg), g = 9,8 m/s2 y c es un coe�ciente que depende del tipode paracaídas y de otros factores, que en este caso supondremos igual a 14. Determinecuál es la masa aproximada del paracaidista si a los 7 segundos de lanzarse del avióncae a una velocidad de 35 m/s.

24. [Mathematica] Dada la ecuación

ln x+2

x2 + 1= 1,

compruebe que el método de Newton-Raphson (2.7) no converge para el valor inicialx0 = 40. Encuentre otro valor inicial para que dicho método converja. Aplique otrométodo que acelere la convergencia a la raíz. Determine la multiplicidad de la raíz.

25. [Mathematica] Aproxime el único cero de la función f(x) = xe−x + 1 mediante elmétodo de Newton-Raphson, con la condición de parada |xn+1 − xn| < 10−4. Compareeste resultado con el obtenido utilizando instrucciones directas deMathematica. ¾Puedeexplicar lo que ocurre si en el método de Newton-Raphson se toma como aproximacióninicial x0 = 4?

26. La función f(x) = x3+x2−5x+3 tiene una raíz doble en x = 1. Aproxime esa soluciónutilizando cuatro iteraciones de los métodos de Newton-Raphson (2.7) y (2.8). Expliquelas diferencias que encuentra al utilizar ambos métodos.

27. [Mathematica] Se sabe que el polinomio p(x) = 100−75x−15x2+30x3−6x4−3x5+x6

tiene dos raíces reales dobles. Aproxime ambas soluciones tanto con el método deNewton-Raphson (2.7) como con el método de Newton-Raphson modi�cado hasta quela diferencia de dos aproximaciones consecutivas sea menor que 10−4. Comente losresultados.

44

28. a) Obtenga las dos primeras aproximaciones (x1, y1) y (x2, y2) de una de las solu-ciones (x, y) del sistema

12(x2 − y2)− 24x+ 7 = 0,

9(x− 1)y + 1 = 0,

}que proporciona el método de Newton-Raphson si se parte del punto (x0, y0) =

(0, 0). Repita el ejercicio para (x0, y0) = (2,−1).b) [Mathematica] Obtenga de forma aproximada mediante el mismo método la solu-

ción de dicho sistema que está en el dominio [0, 2/5] × [0, 2/5], de modo que ladistancia de las dos últimas aproximaciones sea menor que 10−8. Obtenga tambiénuna aproximación de cualquier otra solución real del sistema que se sitúe fuera deese dominio.

29. a) Obtenga la primera aproximación que proporciona el método de Newton-Raphsonde una solución del sistema

x2z − 2y tg xz = 0,

2(x+ 3y)ez−1 + z = 0,

x− 3y − 6z + 4 = 0,

si se parte del punto (x0, y0, z0) = (1,−1,−1).

b) [Mathematica] Aplicando el mismo método, aproxime al menos dos solucionesreales de dicho sistema de modo que la distancia de las dos últimas aproximacionesde cada solución sea menor que 10−6.

30. a) Realizando dos iteraciones del método de Newton-Raphson, obtenga una aproxi-mación de la solución del sistema

x2 sen xy − y = 0,

2y2 ln(x+ 3y) + x = 0,

}situada en las proximidades del punto (−10, 4). Sin realizar más iteraciones, ¾cómose podría medir cuánto se aproxima el resultado obtenido a una solución delsistema?

b) [Mathematica] Aplicando el mismo método, aproxime tres soluciones reales dedicho sistema de modo que la distancia de las dos últimas aproximaciones decada solución sea menor que 10−7.

31. [Mathematica] Observe qué respuesta da la instrucción �NSolve� cuando intenta resolverla ecuación x2 = 1 + senx. Utilice la instrucción �Plot� para saber cuántas solucionestiene la ecuación anterior. Encuentre una aproximación de cada una de esas solucionesutilizando la instrucción �FindRoot�.

32. a) Realizando dos iteraciones del método de Newton-Raphson, obtenga una aproxi-mación de la solución del sistema

x2 − 2x− sen y = −1,

e2y(x− 1) = −1.

}situada en las proximidades del punto (0, 0).

45

b) [Mathematica] Utilice la instrucción �FindRoot� para encontrar algunas solucionesde dicho sistema. Por ejemplo, tome como puntos iniciales los puntos (k, k), conk = 0, 2, 4, 10, 12 (Indicación: en uno de estos valores pueden producirse problemasde convergencia). Pruebe también con otros puntos iniciales (x, y) en los quex = y.

33. Sea f la función de�nida en R mediante f(x) = 3x4 + 4x3 − 102x2 + 180x− 50.

a) Determine el número de raíces reales que posee.

b) Localice la raíz de f con menor valor absoluto en un intervalo de amplitud 1/4.

34. Una función f , cuya grá�ca aparece por duplicado en la �gura 2.2, tiene �como seobserva� un cero en el intervalo [0, 4]. Represente aproximadamente en el grá�co dela izquierda los puntos x1, x2 y x3 de la iteración del método de la regula falsi. En elgrá�co de la derecha, si x−1 = 0 y x0 = 4, haga lo mismo para el método de la secante.

1 2 3 4

-1.5

-1

-0.5

0.5

1 2 3 4

-1.5

-1

-0.5

0.5

Figura 2.2: Aplique los métodos de la regula falsi y de la secante

35. Pruebe que la ecuación x+ arctanx = 1 tiene una única solución real s. Encuentre elintervalo de la forma [z, z + 1] con z ∈ Z que contenga a s. Obtenga las dos primerasaproximaciones de s que proporciona el método de la secante aplicado a ese intervalo.

36. Compruebe, aplicando el teorema de Bolzano, que en el intervalo [−4, 8] hay una raízs del polinomio −72+ 12x− 6x2 + x3. Aplique el método de Newton-Raphson a partirdel punto medio de dicho intervalo. Si hubiera algún problema, resuélvalo de modo quese pueda aplicar este método para aproximar s. Haga entonces alguna iteración.

37. Una persona necesita un préstamo de 150000 euros para la compra de un inmueble,préstamo que desea amortizar mediante pagos que no superen los 12000 euros anuales,durante 25 años. La fórmula de amortización del préstamo es

150000 =100M

i

[1−

(100

100 + i

)25],

donde M es la cantidad a pagar cada año e i es el interés anual del préstamo (en%).¾Cuál es el interés anual máximo imax que puede asumir?

a) Localice en un intervalo de amplitud 4 unidades el valor de imax.

46

b) Obtenga la fórmula del método de Newton-Raphson para resolver este problema,elija un valor inicial adecuado y obtenga una primera aproximación de imax.

38. Obtenga una aproximación de una solución del sistema

x2 + ey = 5,

sen x+ cos y = 1,

}partiendo del punto (0, π/2) como primera aproximación de la solución, y obteniendosolo la primera iteración que proporciona el método de Newton-Raphson para sistemas.Sin realizar más iteraciones, ¾cómo se podría medir cuánto se aproxima el resultadoobtenido a una solución del sistema?

Algunas cuestiones teóricas

1. De�na el concepto de raíz de multiplicidad k de una ecuación.

2. Si hubiese que aproximar la raíz de multiplicidad 4 del ejercicio 1 por alguno de losmétodos de Newton-Raphson, ¾cuál elegiría? ¾Por qué?

3. Explique qué métodos aplicaría y por qué, para aproximar con el menor número deiteraciones posible los dos ceros de una función f cuya grá�ca es la representada enla �gura 2.3. Exprese la fórmula para la primera iteración en cada caso, eligiendo lospuntos de partida que estime más convenientes.

0.5 1 1.5

-0.2

-0.1

0.1

0.2

Figura 2.3: ¾Qué métodos emplearía para aproximar los ceros de esta grá�ca?

4. La grá�ca de una función aparece por duplicado en la �gura 2.4. Represente aproxi-madamente en el grá�co de la izquierda los primeros puntos de la iteración del métodode la regula falsi, partiendo de los extremos del intervalo donde se ha representado lagrá�ca. En el grá�co de la derecha haga lo mismo para el método de la secante. ¾Haciaqué solución se tiende en cada caso?

5. ¾Qué diferencias observa entre los métodos de Newton-Raphson, por un lado, y losmétodos de bisección, de la regula falsi y de la secante, por otro?

6. Determine la verdad o falsedad de las siguientes a�rmaciones.

47

0.5 1 1.5

-0.2

-0.1

0.1

0.2

0.5 1 1.5

-0.2

-0.1

0.1

0.2

Figura 2.4: Aplique los métodos de la regula falsi y de la secante

a) Si α es una raíz triple de la ecuación f(x) = 0, entonces α es una raíz doble de laecuación f ′(x) = 0.

b) El método de la secante no siempre converge.

c) El método de la regula falsi no siempre converge.

d) El método de la secante puede utilizarse para aproximar raíces múltiples (aunquequizá lo haga con lentitud).

e) Si p(x) es un polinomio de grado k, su derivada k-ésima p(k)(x) es una funciónconstante. Por tanto, p(x) no puede tener raíces de multiplicidad k + 1.

7. Represente una grá�ca que corte al eje x en un punto, y aplique sobre la grá�ca elmétodo de la secante para aproximar el cero de la función. Repita el ejercicio (engrá�ca aparte similar) para el método de Newton-Raphson.

8. Comente las ventajas e inconvenientes que tiene el método de la secante respecto delmétodo de bisección.

48

Capítulo 3

RESOLUCIÓN NUMÉRICADE SISTEMAS DE ECUACIONESLINEALES

En este capítulo indagaremos sobre diversas cuestiones acerca de la resolución de sistemaslineales compatibles determinados. Especialmente nos centraremos en su resolución numérica,que es clave para la resolución de sistemas de gran tamaño. También introduciremos unmétodo para la obtención aproximada de valores propios de una matriz.

3.1. Preliminares de Álgebra Lineal

En esta sección recordamos algunos conceptos y resultados de la asignatura de Álge-bra Lineal, que necesitaremos a lo largo de este capítulo; e introducimos otros conceptos yresultados de esa rama de las Matemáticas que tendrán relevancia en las siguientes secciones.

Representaremos por Mm×n(R) al conjunto de las matrices reales m × n (m �las y n

columnas). Del mismo modo, Mm×n(C) representa a las matrices complejas.1 A los ele-mentos de estos conjuntos los representaremos frecuentemente por (aij) (i = 1, 2, . . . ,m,j = 1, 2, . . . , n): si escribimos A = (aij) ∈ Mm×n(R) querrá decir que

A =

a11 a12 · · · a1na21 a22 · · · a2n...

......

am1 am2 · · · amn

, (3.1)

donde aij ∈ R para cualesquiera i = 1, 2, . . . ,m y j = 1, 2, . . . , n. Si m = 1, entonces se diceque (aij) es una matriz �la; si n = 1, se dice que (aij) es una matriz columna; si m = n,se dice que (aij) es una matriz cuadrada de orden n. Al conjunto formado por estas últimasse le representará más simpli�cadamente por Mn(R) (si son reales) o por Mn(C) (si soncomplejas).

1En lo que sigue, cuando introduzcamos conceptos y resultados sobre matrices, entenderemos que las

matrices pueden ser reales o complejas salvo que se especi�que otra cosa. En ocasiones hablaremos de

matrices reales no porque ese concepto o resultado no pueda extenderse a matrices complejas, sino porque

solo lo necesitaremos para matrices reales.

49

Podemos apreciar que Rn es un conjunto bastante similar a M1×n(R), y también aMn×1(R) (algo semejante puede decirse de Cn). De hecho, estos tres conjuntos son lo que enálgebra se denomina espacios isomorfos : sus elementos y sus operaciones son esencialmentelos mismos. Esto permite identi�carlos y usarlos casi indistintamente. Por ejemplo, muchasveces interesará representar los elementos de Rn como matrices columna n × 1 o 1 × n. Dehecho, frecuentemente nos referiremos a estas matrices como puntos o vectores de Rn, yviceversa.

Los elementos de una matriz �la o de una matriz columna se denotan usualmente con unsolo subíndice, puesto que el otro se sobrentiende que es igual a 1: por ejemplo,

a = (a1 a2 · · · an) , b =

b1b2...bn

.

Sea A = (aij) una matrizm×n. Se de�ne entonces su matriz traspuesta, que denotaremospor At, como la matriz que se obtiene a partir de A intercambiando sus �las por sus columnasen el mismo orden, esto es,

At =

a11 a12 · · · a1n...

......

am1 am2 · · · amn

t

=

a11 · · · am1

a12 · · · am2

......

a1n · · · amn

.

Una matriz cuadrada A se dice que es simétrica si coincide con su traspuesta: At = A.Dado un sistema lineal cualquiera

a11x1 + a12x2 + · · ·+ a1nxn = b1,

a21x1 + a22x2 + · · ·+ a2nxn = b2,

· · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·am1x1 + am2x2 + · · ·+ amnxn = bm,

(3.2)

la matriz de coe�cientes del sistema es la matriz A dada por (3.1) y la matriz de los términos

independientes del sistema es la matriz columna b = (b1 b2 . . . bm)t. Si x = (x1 x2 . . . xn)

t esla matriz columna formada por las incógnitas del sistema, entonces el sistema (3.2) puedeexpresarse de la forma Ax = b, esto es:

a11 a12 · · · a1na21 a22 · · · a2n...

......

am1 am2 · · · amn

x1

x2

...xn

=

b1b2...bm

.

Si A = (aij) es una matriz cuadrada de orden n, su diagonal principal es la formada porlos elementos aii, con i = 1, 2, . . . , n. Una matriz cuadrada A = (aij) se dice que es triangularsuperior si los elementos que se sitúan bajo la diagonal principal son todos nulos, esto es, siaij = 0 para todo i > j. Análogamente, A es triangular inferior si aij = 0 para todo i < j. Si

50

una matriz A de orden n es triangular superior e inferior a la vez �es decir, si aij = 0 paratodo i = j�, entonces se dice que A es una matriz diagonal. En tal caso la matriz puededenotarse de la forma

A = diag (a11, a22, . . . , ann) =

a11 0 · · · 0

0 a22 · · · 0...

.... . .

...0 0 · · · ann

.

El ejemplo más característico de matriz diagonal de orden n, es la matriz identidad, quedenotaremos por In, y que es la matriz en la que todos los elementos de su diagonal principalson iguales a 1 (y el resto son cero):

In = diag (1, 1, n. . ., 1) =

1 0 · · · 0

0 1 · · · 0...

.... . .

...0 0 · · · 1

.

Una matriz cuadrada A = (aij) de orden n se dice que es diagonalmente dominante por

�las en sentido estricto si el módulo de cada elemento de la diagonal principal es mayor quela suma de los módulos de los demás elementos de su �la, es decir, si

|aii| >n∑

j=1j =i

|aij| para todo i = 1, 2, . . . , n.

Del mismo modo, se dice que A es diagonalmente dominante por columnas en sentido estricto

si el módulo de cada elemento de la diagonal principal es mayor que la suma de los módulosde los demás elementos de su columna, esto es, si

|ajj| >n∑

i=1i =j

|aij| para todo j = 1, 2, . . . , n.

A veces, por brevedad, omitiremos las palabras �en sentido estricto� al referirnos a ambasde�niciones, ya que en esta asignatura no se utilizan las correspondientes de�niciones noestrictas.

Ejemplo 3.1 Dadas las matrices

A =

−9 −1 2 1

−1 5 −3 0

1 1 −4 1

1 −2 1 5

, B =

4 −1 2 1

1 −5 −3 0

1 1 7 1

1 −2 1 5

y C =

5 −1 2 1

1 −5 π 0

2 −1 7 1

1 −e 1 9

,

es claro que A es diagonalmente dominante por �las en sentido estricto pero no por columnas(en su tercera columna ocurre que | − 4| < |2|+ | − 3|+ |1|), B es diagonalmente dominantepor columnas en sentido estricto pero no por �las (en su primera �la ocurre que |4| =

| − 1| + |2| + |1|) y C es diagonalmente dominante tanto por �las como por columnas ensentido estricto. 2

51

Al determinante de una matriz cuadrada A lo denotaremos |A| o det(A). Una matrizcuadrada A de orden n es inversible cuando existe otra A−1 tal que AA−1 = A−1A = In.Equivalentemente, A es inversible si y solo si |A| = 0. El siguiente resultado muestra unprimer tipo de matrices que siempre son inversibles.

Teorema 3.1 Toda matriz diagonalmente dominante (por �las o por columnas) en sentido

estricto, es inversible.

Sea A una matriz cuadrada real. Se dice que λ ∈ C es un valor propio de la matriz A

si existe un vector x ∈ Cn no nulo tal que Ax = λx. En este caso se dice también que x

es un vector propio de A asociado al valor propio λ. Entre los vectores propios asociados acualquier valor propio λ se incluye también al vector θ = (0, 0, . . . , 0), puesto que Aθ = λθ.

Si λ es un valor propio de una matriz, al conjunto formado por todos los vectores propiosasociados a λ lo denotaremos Vλ. Sin embargo, en el caso particular de que λ sea un valorpropio real, solo utilizaremos los vectores de Vλ cuyas componentes sean reales, es decir, losque pertenecen a Rn; al conjunto formado por estos vectores propios lo denotaremos Wλ.2

El siguiente teorema muestra cómo pueden obtenerse los valores propios de una matrizcuadrada así como los vectores propios asociados a cada valor propio.

Teorema 3.2 Para cualquier matriz cuadrada real A de orden n se satisfacen las siguientes

propiedades:

1. λ ∈ C es un valor propio de A si y solo si |A− λIn| = 0.

2. Si λ es un valor propio de A, entonces el conjunto Vλ de sus vectores propios coincide

con el conjunto de las soluciones complejas del sistema (A− λIn)x = 0. En el caso

de que λ sea real, Wλ se obtiene como el conjunto de las soluciones reales de dicho

sistema.

3. Para cada valor propio λ de A, se tiene que Vλ es un subespacio vectorial de Cn.

En el caso de que λ ∈ R, entonces Wλ es un subespacio vectorial de Rn; además

dimWλ = dimVλ, de modo que cualquier base de Wλ es también una base de Vλ.

Para la aplicación del primer apartado del teorema 3.2, observe que el determinante|A− λIn| es un polinomio respecto de la variable λ, p(λ), al que se denomina polinomio

característico de la matriz A. Es claro que p(λ) es un polinomio con coe�cientes reales, y quesu grado es n. Luego los valores propios de A son las raíces de su polinomio característico.A la multiplicidad rλ de cada raíz λ se le denomina multiplicidad algebraica del valor propioλ.3

Otra clara consecuencia del teorema 3.2 es que, para cada valor propio λ de una matrizA de orden n, el sistema homogéneo (A− λIn)x = 0 es compatible indeterminado, de dondese deduce que la dimensión sλ de Vλ (y la de Wλ, en el caso λ ∈ R) es mayor o igual que 1. Asλ se le denomina multiplicidad geométrica de λ. Como consecuencia de resultados conocidos

2En resumidas cuentas, Vλ ⊂ Cn para todo valor propio λ, real o imaginario. Y si λ ∈ R, entoncesWλ ⊂ Rn y puede probarse que Wλ Vλ: muchos vectores de Vλ no pertenecen a Rn.

3Por tanto, la suma de las multiplicidades algebraicas de los valores propios de una matriz de orden n es

siempre igual a n.

52

del Álgebra lineal se tiene que sλ = n− rg(A−λIn). Además, el siguiente resultado establecela relación que existe entre las multiplicidades algebraica y geométrica.

Teorema 3.3 Sea A una matriz cuadrada real de orden n, y sea λ cualquiera de sus valores

propios. Entonces se cumple que 1 ≤ sλ ≤ rλ.

El radio espectral de una matriz cuadrada A de orden n se de�ne como el máximo de losmódulos de sus valores propios; se denota por ρ(A). Todo valor propio λ de A cuyo módulosea máximo, es decir, que veri�que que |λ| = ρ(A), se dice que es dominante. Si hay un únicovalor propio dominante (teniendo en cuenta también su multiplicidad, que deberá ser 1), sedice que es estrictamente dominante.

En el caso de matrices reales, el valor propio estrictamente dominante, si existe, tieneque ser real (si fuera imaginario, su conjugado sería otro valor propio con el mismo módulo,por lo que no sería estrictamente dominante).

Es fácil probar que el polinomio característico de una matriz cuadrada A coincide conel de su traspuesta. Por tanto A y At tienen los mismos valores propios y, en particular,ρ(A) = ρ(At).

Sea A una matriz cuadrada real de orden n. Se dice que A es diagonalizable cuandoes semejante a una matriz diagonal, es decir, cuando existe una matriz P inversible y unamatriz D diagonal tales que P−1AP = D (las matrices P y D podrán ser complejas).

El siguiente teorema proporciona una caracterización del concepto de matriz diagonali-zable.

Teorema 3.4 Sea A una matriz cuadrada real de orden n.

1. A es diagonalizable si y solo si las multiplicidades algebraica y geométrica de cada uno

de sus valores propios coinciden: rλ = sλ para todo λ valor propio de A (usando la

notación del teorema 3.3).

2. Si A es diagonalizable y λ1, λ2, . . . , λn son sus valores propios escritos en cualquier

orden y apareciendo repetidos tantas veces como indique su multiplicidad, entonces A

es semejante a la matriz D = diag(λ1, λ2, . . . , λn). Además, si para cada i = 1, 2, . . . , n

se elige un vector propio xi asociado a λi de manera que los vectores propios elegidos

asociados a λi formen una base de su subespacio la matriz P tal que P−1AP = D

siempre que los vectores columna x1, x2, . . . , xn de P sean n vectores propios linealmente

independientes de A respectivamente asociados a los valores propios λ1, λ2, . . . , λn.

Como consecuencia inmediata del teorema anterior, se obtiene el siguiente resultado.

Teorema 3.5 Si una matriz cuadrada real de orden n posee n valores propios distintos,

entonces es diagonalizable.

Sea A una matriz cuadrada de orden n. Una matriz se dice que es la submatriz principal

de orden k de A, con k = 1, 2, . . . , n, si se obtiene a partir de A eliminando sus últimas n−k

�las y n−k columnas. A esta submatriz la denotaremos A(k). Observe que A(n) = A. Se de�neel menor principal de orden k de la matriz A, denotado por mk(A), como el determinantede la submatriz principal de orden k de la matriz A: mk(A) = |A(k)|.

53

Sea A una matriz simétrica real de orden n. Se dice que A es de�nida positiva si todossus menores principales son positivos: mk(A) > 0 para todo k = 1, 2, . . . , n. En particular, eldeterminante de una matriz de�nida positiva es positivo, y por tanto la matriz es inversible.

El siguiente teorema proporciona dos propiedades de las matrices simétricas reales y dosmás de las que son de�nidas positivas.

Teorema 3.6 Para cualquier matriz simétrica real A se cumplen las siguientes propiedades:

1. Todos los valores propios de A son reales.

2. A es diagonalizable.

3. A es de�nida positiva si y solo si todos sus valores propios son positivos.

4. Si A es de�nida positiva, entonces todos los elementos de su diagonal son positivos.

3.2. Normas matriciales

En la sección 1.3 se introdujo el concepto de espacio normado, y se comentó que unprimer ejemplo de espacio normado era el espacio euclídeo n-dimensional. En esta secciónse introducen otros espacios normados y, en particular, las normas matriciales y algunosconceptos relacionados con ellas. Esta sección, como la anterior, proporciona herramientasbásicas para los métodos numéricos que se introducirán en este capítulo.

En el conjunto Rn, además de la norma euclídea pueden de�nirse un número in�nito denormas. En concreto, se puede probar que

∥u∥p =

(n∑

i=1

|ui|p)1/p

de�ne una norma en Rn para todo p ∈ [1,∞). Observe que la norma euclídea se obtienecuando p = 2. También nos interesará de modo especial el caso p = 1:

∥u∥1 =n∑

i=1

|ui|.

Otra norma que utilizaremos en Rn es la de�nida como sigue:

∥u∥∞ = max1≤i≤n

|ui|.

Por contraste con las normas que de�niremos a continuación, a las normas de�nidas enRn las denominaremos normas vectoriales.

En el conjunto Mn(R) pueden de�nirse muchas normas diferentes que doten a dichoconjunto de estructura de espacio normado. Pero solo nos interesará trabajar con normasque cumplan las propiedades que introduciremos en las de�niciones 3.1 y 3.3.

De�nición 3.1 Una norma ∥ · ∥ en Mn(R) se dice que es una norma matricial si cumplela siguiente propiedad añadida:

∥AB∥ ≤ ∥A∥∥B∥ para cualesquiera A,B ∈ Mn(R). (3.3)

54

Observe que una norma matricial cumple la siguiente propiedad:

∥In∥ ≥ 1. (3.4)

La razón es simple: aplicando (3.3), se tiene que ∥In∥ = ∥InIn∥ ≤ ∥In∥∥In∥ = ∥In∥2, dedonde se deduce trivialmente la propiedad.

A continuación introducimos las tres normas matriciales más utilizadas.

De�nición 3.2 Podemos introducir tres normas matriciales en el conjunto Mn(R) de lasiguiente forma. Para cada A = (aij) ∈ Mn(R) se de�ne:

∥A∥2 =√ρ(AtA),

∥A∥1 = max1≤j≤n

n∑i=1

|aij|,

∥A∥∞ = max1≤i≤n

n∑j=1

|aij|.

Ejemplo 3.2 Es inmediato que

∥In∥1 = ∥In∥2 = ∥In∥∞ = 1.

Por ejemplo, como el único valor propio de In es λ = 1, entonces ∥In∥2 =√ρ(I tnIn) =√

ρ(In) =√1 = 1. Las otras dos igualdades son triviales. Por tanto, en este triple ejemplo

se tiene la igualdad en (3.3). Pero hay también muchas normas matriciales para las que∥In∥ > 1. 2

El siguiente teorema muestra algunas propiedades del radio espectral y de las normasmatriciales.

Teorema 3.7 Sea A una matriz cuadrada real. Entonces:

1. El radio espectral de A es menor o igual que cualquiera de sus normas matriciales:

ρ(A) ≤ ∥A∥M para cualquier norma matricial ∥ · ∥M . Es más, el radio espectral de A

es el ín�mo de todas las normas matriciales de A:

ρ(A) = ınf {∥A∥M : ∥ · ∥M es una norma matricial} .

2. En particular, ρ(A) < 1 si y solo si existe una norma matricial ∥·∥M tal que ∥A∥M < 1.

De�nición 3.3 Sea ∥ · ∥M una norma matricial en Mn(R), y sea ∥ · ∥V una norma vectorialen Rn. Se dice que estas normas son compatibles si para cada matriz A ∈ Mn(R) y cadavector u ∈ Rn (expresado como matriz columna) se tiene que

∥Au∥V ≤ ∥A∥M∥u∥V .

Puede probarse que las normas matriciales ∥ · ∥1, ∥ · ∥2 y ∥ · ∥∞ son respectivamentecompatibles con las normas vectoriales denotadas del mismo modo.

55

3.3. Condicionamiento de una matriz

En esta sección se introduce un concepto, el de condicionamiento de una matriz realinversible, que es importante tener en cuenta para la resolución de sistemas lineales quepuedan tener errores (incluso muy pequeños) en sus coe�cientes o términos independientes.Por brevedad solo consideraremos el caso de errores en los términos independientes.

De�nición 3.4 Fijada una norma matricial ∥ · ∥ en Mn(R), para cada matriz inversibleA ∈ Mn(R) se de�ne su condicionamiento mediante la siguiente expresión:

cond(A) = ∥A∥∥A−1∥.

Por tanto, para cada matriz real inversible A, nosotros podremos considerar tres de�ni-ciones de condicionamiento, una para cada una de las normas matriciales introducidas:

cond1(A) = ∥A∥1∥A−1∥1, cond2(A) = ∥A∥2∥A−1∥2, cond∞(A) = ∥A∥∞∥A−1∥∞.

El siguiente teorema muestra una propiedad del condicionamiento de una matriz.

Teorema 3.8 Sea A una matriz real inversible de orden n, y sea ∥ · ∥ cualquier norma

matricial en Mn(R). Entonces se tiene que

cond(A) ≥ ∥In∥ ≥ 1.

El condicionamiento de una matriz no depende de la magnitud de sus coe�cientes: porejemplo, dada una matriz real inversible A, es fácil probar que cond(kA) = cond(A) paratodo k = 0.

Ejemplo 3.3 Como ∥In∥1 = ∥In∥2 = ∥In∥∞ = 1 (recuerde el ejemplo 3.2), es claro que

cond1(In) = cond2(In) = cond∞(In) = 1.

En general, para cualquier matriz diagonal A = diag(a11, a22, . . . , ann) inversible (esto es, talque aii = 0 para todo i = 1, 2, . . . , n), se puede probar que

cond1(A) = cond2(A) = cond∞(A) =max1≤i≤n

|aii|

mın1≤i≤n

|aii|.

El lector podrá, al menos, comprobar este último resultado para algún ejemplo concreto. 2

A continuación trataremos acerca de los efectos que puede tener el condicionamiento de lamatriz de coe�cientes en la resolución de un sistema compatible determinado. Mostraremosque el condicionamiento de dicha matriz puede afectar a la exactitud de la resolución delsistema. Pero nos limitaremos a estudiar un supuesto sencillo, el que muestra el siguienteteorema.

56

Teorema 3.9 Sea A una matriz cuadrada real inversible de orden n, y sea b una matriz

columna real del mismo orden. Supongamos que de la matriz b solo se conoce una aproxi-

mación b∗. Sea x (= A−1b) la solución del sistema Ax = b y x∗ (= A−1b∗) la solución del

sistema Ax = b∗. Elegimos una norma matricial en Mn(R) y una norma vectorial en Rn

que sean compatibles. Si consideramos los errores relativos

rb =∥b− b∗∥

∥b∥y rx =

∥x− x∗∥∥x∥

,

entonces se tiene la siguiente desigualdad:rb

cond(A)≤ rx ≤ rb cond(A). (3.5)

En conclusión, la segunda desigualdad en (3.5) puede ser problemática si el condi-cionamiento de la matriz A es un número grande. Teóricamente, el error relativo en b podríamultiplicarse en la solución x por el valor del condicionamiento.

Ejemplo 3.4 Consideramos el sistema compatible determinado

−x+ 2y + z = 2,

150x− 300y − 151z = −300,

250x− 500,1y − 250,9z = −500,

en el que los términos independientes anteriores son una aproximación de los verdaderos.Siguiendo la notación utilizada, tenemos pues que

A =

−1 2 1

150 −300 −151

250 −500,1 −250,9

, b∗ =

2

−300

−500

,

y la solución del sistema es x∗ = A−1b∗ = (−2 0 0)t.Es fácil calcular que el condicionamiento de la matriz A con respecto a las normas ∥ · ∥1

y ∥ · ∥∞ es aproximadamente igual a 3 · 106 y 2,49 · 106, respectivamente.Si los términos independientes exactos son los de la matriz b = (2 −301 −501)t, entonces

la solución verdadera es x = A−1b = (1 1 1)t. Se observa que un error no muy grande en lostérminos independientes ha cambiado sustancialmente la solución del sistema.

Utilizando las normas del in�nito, el error relativo en los términos independientes es

rb =∥b− b∗∥∞

∥b∥∞=

∥(0, 1, 1)∥∞∥(2,−301,−501)∥∞

=1

501,

mientras que el error relativo en la solución es

rx =∥x− x∗∥∞

∥x∥∞=

∥(3, 1, 1)∥∞∥(1, 1, 1)∥∞

= 3.

Luego rx = 1503 · rb: el error relativo en la solución ha multiplicado 1503 veces el errorrelativo en los términos independientes (y podría haber sido peor: el condicionamiento deA indica que el factor multiplicativo podría haber sido de hasta 2,49 · 106). De hecho, sihubiéramos tomado b∗ = (2,01 −301,01 −500,99)t �como se pide en el ejercicio 4 de lasección 3.7�, entonces el resultado es aún más sorprendente, puesto que se obtiene

x∗ = (−23,88 −10,69 −0,49)t, rb ≃ 1,99601 · 10−5 y rx ≃ 24,88,

con lo que rx resulta ser más de 1,2 millones de veces mayor que rb. 2

57

3.4. Métodos iterativos de resolución de sistemas lineales

Con esta sección comenzamos el estudio del principal objetivo de este capítulo: encontrarmétodos que puedan proporcionar una aproximación precisa de la solución de un sistemalineal n×n compatible determinado. Estos métodos se podrán aplicar a sistemas lineales decualquier orden, pero su interés reside en su aplicación a sistemas de un orden tan elevadoque su resolución mediante métodos exactos pueda producir excesivos errores de redondeoen la solución. Por el contrario, los errores de redondeo pierden su in�uencia en los métodositerativos por el modo en que estos están diseñados.

Así, tanto en esta sección como en la siguiente trabajaremos con sistemas lineales Ax = b,donde A es una matriz cuadrada real inversible de orden n (y b ∈ Rn). Los métodos iterativosde resolución de estos sistemas proporcionan una sucesión de puntos de Rn que convergen asu solución. Aunque no lo demostremos, estos métodos iterativos siempre pueden expresarsemediante una fórmula matricial del siguiente tipo:

x(k+1) = Cx(k) + d, k = 0, 1, 2, . . . , (3.6)

donde C ∈ Mn(R) y d ∈ Rn se de�nen a partir de las matrices A y b del sistema, y el puntode partida x(0) puede ser cualquier punto de Rn.4 El siguiente resultado precisa cuándo lasiteraciones del tipo (3.6) convergen a un punto.5

Teorema 3.10 Si C es una matriz cuadrada real de orden n, y d y x(0) son dos puntos de

Rn, la iteración (3.6) converge si y solo si ρ(C) < 1.

Por tanto, que una iteración del tipo (3.6) converja solo depende de la matriz C: si suradio espectral es menor que 1, los métodos de la forma (3.6) convergen para cualesquiera d

y x(0) en Rn.El siguiente resultado podrá aplicarse para acotar el error de aproximación de los métodos

iterativos de resolución de sistemas.

Teorema 3.11 Sea C una matriz cuadrada real de orden n tal que ρ(C) < 1, y sean d y

x(0) dos puntos de Rn. Si {x(k)} es la sucesión de puntos dada por (3.6) y x = lımk→∞

x(k),6

entonces

∥x(k) − x∥V ≤ ∥C∥M1− ∥C∥M

∥x(k) − x(k−1)∥V , k = 1, 2, . . . (3.7)

para cualquier norma matricial ∥ · ∥M tal que ∥C∥ < 1,7 y para cualquier norma vectorial

∥ · ∥V compatible con ∥ · ∥M .

En las aplicaciones de este teorema resultará que el término de la izquierda en (3.7) esdesconocido (puesto que no se conoce la solución x), de manera que el término de la derechaproporcionará una cota del error ∥x(k)−x∥V , que es proporcional a la norma de la diferencia

4Como regla general, si el método es convergente, cuanto más próximo esté x(0) a la solución del sistema,

se necesitarán menos iteraciones para llegar a una buena aproximación de la solución.5Nosotros lo aplicaremos al caso en el que ese punto es la solución de un sistema lineal dado, pero el

teorema es válido también cuando C y d no provienen de un sistema.6El teorema 3.10 asegura que existe este límite.7Como ρ(C) < 1, el teorema 3.7 asegura la existencia de normas matriciales ∥ · ∥M tales que ∥C∥M < 1.

58

de las dos últimas aproximaciones obtenidas. Observe también que es indispensable utilizaruna norma matricial tal que ∥C∥M < 1 (si ∥C∥M ≥ 1, el término de la derecha nunca podríaser una cota del error ∥x(k) − x∥V , ya que dicho término sería negativo o in�nito).

3.5. Métodos de Jacobi y de Gauss-Seidel

En esta sección introducimos dos métodos iterativos para la resolución aproximada desistemas lineales Ax = b del tipo establecido en la sección 3.4, y estudiamos en qué casosson aplicables. Ambos métodos (y otros muchos) se basan en la siguiente descomposición dela matriz de coe�cientes A = (aij) del sistema:

A = D − L− U, (3.8)

donde D = diag(a11, a22, . . . , ann), y L = (lij) y U = (uij) son las matrices triangularesde�nidas por

lij =

{0 i ≤ j,

−aij i > j,y uij =

{−aij i < j,

0 i ≥ j.

Por ejemplo, la matriz L tendría el siguiente aspecto:

L =

0 0 · · · · · · · · · 0

−a21 0 0 · · · · · · 0

−a31 −a32 0 0 · · · 0...

. . . . . . . . . . . ....

.... . . . . . . . . 0

−an1 −an2 · · · −ann−2 −ann−1 0

.

Ejemplo 3.5 Si A =

9 −1 2 1

1 2 −3 0

−1 1 0 −1

1 −2 1 3

, su descomposición de la forma (3.8) es

A =

9 0 0 0

0 2 0 0

0 0 0 0

0 0 0 3

−

0 0 0 0

−1 0 0 0

1 −1 0 0

−1 2 −1 0

−

0 1 −2 −1

0 0 3 0

0 0 0 1

0 0 0 0

. 2

Para que se puedan aplicar los métodos de Jacobi y de Gauss-Seidel, será necesario exigirque todos los elementos de la diagonal de la matriz de coe�cientes del sistema sean no nulos.Por ejemplo, si Ax = b es un sistema donde A = (aij) es la matriz del ejemplo 3.5, entonces,como det(A) = 50 = 0, el sistema es compatible determinado. Pero a este sistema no se lepodrá aplicar ni el método de Jacobi ni el de Gauss-Seidel, ya que a33 = 0. Obviamente,si se intercambia la posición de las ecuaciones tercera y cuarta, se obtendría un sistemaequivalente A′x = b′ en el que a′ii = 0 para todo i = 1, 2, 3, 4; con lo que se podrá intentaraplicar ambos métodos para este último sistema.

59

En general, los métodos de Jacobi y de Gauss-Seidel son sensibles al cambio del orden delas ecuaciones: puede que converjan para cierto orden determinado a la solución del sistema,y que para otro orden no converjan. Ahondaremos en este punto a continuación, al introducirambos métodos, y también mediante los ejercicios propuestos.

3.5.1. Método de Jacobi

Sea Ax = b un sistema del tipo establecido en la sección 3.4, y sea A = D − L − U ladescomposición de A introducida en esa sección. Entonces, el método de Jacobi para dichosistema se de�ne mediante la iteración

x(k+1) = D−1(L+ U)x(k) +D−1b, k = 0, 1, 2, . . . , (3.9)

que parte de cualquier punto x(0) ∈ Rn. Es evidente que la iteración de Jacobi es del tipo(3.6), donde C = D−1(L+U) y d = D−1b. A la primera de estas matrices la denominaremoscomo la matriz del método de Jacobi y, en cuanto que solo depende de la matriz A (y de ladescomposición (3.8) realizada en ella), la denotaremos por JA:

JA = D−1(L+ U). (3.10)

Observe que JA = D−1(L + U) = D−1(D − A) = D−1D − D−1A = In − D−1A. Así,JA = In − D−1A es una de�nición equivalente de la matriz del método de Jacobi. En losejercicios, usar una u otra expresión de JA es indiferente. Desarrollando cualquiera de esasexpresiones de JA y realizando también el producto D−1b, se obtiene

JA =

0−a12a11

−a13a11

· · · −a1na11

−a21a22

0−a23a22

· · · −a2na22

.... . . . . . . . .

...

−an−11an−1n−1

−an−12an−1n−1

. . . 0−an−1nan−1n−1

−an1ann

−an2ann

· · · −ann−1

ann0

, D−1b =

b1a11b2a22...

bnann

.

Sustituyendo estos resultados en (3.9) e igualando componente a componente, se obtiene unafórmula equivalente de la iteración del método de Jacobi:

x(k+1)i =

1

aii

(bi −

i−1∑j=1

aijx(k)j −

n∑j=i+1

aijx(k)j

), i = 1, 2, . . . , n, k = 0, 1, 2, . . . , (3.11)

donde los sumatorios con contenido nulo �los que van desde j = 1 hasta j = 0, o bien desdej = n+ 1 hasta j = n� se entiende que son iguales a cero.

El siguiente teorema determina cuándo converge el método de Jacobi.

60

Teorema 3.12 Sea A una matriz real inversible de orden n tal que aii = 0 para todo i =

1, 2, . . . , n; y sea b una matriz columna real del mismo orden. Consideramos la descomposi-

ción de A dada por (3.8) y la matriz JA de�nida por (3.10). Y para cualquier x(0) ∈ Rn, sea

{x(k)} la sucesión de puntos de Rn obtenida mediante la iteración (3.9) o (3.11). Entonces{x(k)} converge a la solución del sistema Ax = b si y solo si ρ(JA) < 1.

Los dos siguientes resultados proporcionan condiciones su�cientes (aunque no necesarias,a diferencia del teorema 3.12) para que converja el método de Jacobi. Son muy útiles enbastantes ocasiones.

Teorema 3.13 Sea A una matriz real de orden n, y sea b una matriz columna real del

mismo orden. Si la matriz A es diagonalmente dominante (por �las o por columnas) en

sentido estricto,8 entonces la iteración de Jacobi (3.9) o (3.11) converge a la solución del

sistema Ax = b para cualquier punto de partida x(0) ∈ Rn.

Teorema 3.14 Sea A una matriz simétrica real de orden n, y sea b una matriz columna

real del mismo orden. Consideramos la descomposición de A dada por (3.8). Si las matrices

A(= D−L−U) y D+L+U son de�nidas positivas, entonces la iteración de Jacobi (3.9) o(3.11) converge a la solución del sistema Ax = b para cualquier punto de partida x(0) ∈ Rn.9

3.5.2. Método de Gauss-Seidel

Sea Ax = b un sistema del tipo establecido en la sección 3.4, y sea A = D − L − U ladescomposición de A introducida en esta sección. Entonces, el método de Gauss-Seidel paradicho sistema se de�ne mediante la iteración

x(k+1) = (D − L)−1Ux(k) + (D − L)−1b, k = 0, 1, 2, . . . , (3.12)

que parte de cualquier punto x(0) ∈ Rn. Luego es una iteración del tipo (3.6), donde C =

(D − L)−1U y d = (D − L)−1b. A la primera de estas matrices se la denomina como lamatriz del método de Gauss-Seidel ; como su de�nición solo depende de la matriz A y de sudescomposición (3.8), la denotaremos por GA:

GA = (D − L)−1U. (3.13)

Observe queGA = (D−L)−1U = (D−L)−1(D−L−A) = (D−L)−1(D−L)−(D−L)−1A =

In − (D−L)−1A. Así, GA = In − (D−L)−1A es una de�nición equivalente de la matriz delmétodo de Gauss-Seidel. Usar una u otra expresión de GA es indiferente.

Observe que la expresión (3.12) es equivalente a la siguiente:

(D − L)x(k+1) = Ux(k) + b, k = 0, 1, 2, . . . , (3.14)

8Observe que, al ser A diagonalmente dominante en sentido estricto, se cumple que aii = 0 para todo

i = 1, 2, . . . , n. Y el teorema 3.1 asegura además que A es inversible.9Recuerde que en la sección 3.1 se vio que el determinante de toda matriz de�nida positiva es positivo;

por tanto, el sistema Ax = b es compatible determinado. Por otra parte, el cuarto apartado del teorema 3.6

asegura que los elementos de la diagonal de A son positivos, luego no nulos.

61

que es un sistema lineal triangular (la matriz de coe�cientes D − L es triangular inferior)con incógnitas x

(k+1)1 , x

(k+1)2 , . . . , x

(k+1)n . Realizando las operaciones matriciales en (3.14) se

obtienei∑

j=1

aijx(k+1)j = bi −

n∑j=i+1

aijx(k)j , i = 1, 2, . . . , n, k = 0, 1, 2, . . . (3.15)

Para cualquier k = 0, 1, 2, . . ., la solución del sistema (3.15) (su resolución es trivial, por serun sistema triangular) es

x(k+1)i =

1

aii

(bi −

i−1∑j=1

aijx(k+1)j −

n∑j=i+1

aijx(k)j

), i = 1, 2, . . . , n, k = 0, 1, 2, . . . (3.16)

Observe que la expresión (3.16) es muy similar a la correspondiente expresión (3.11)del método de Jacobi: la única diferencia está en los superíndices del primer sumatorio (kpara el método de Jacobi y k + 1 para Gauss-Seidel). Así, en el método de Jacobi todaslas componentes de la aproximación x(k+1) se obtienen en función de las componentes dela aproximación anterior x(k). Sin embargo, en el método de Gauss-Seidel, a medida quese van calculando componentes de x(k+1), estas se van utilizando para el cálculo de lassiguientes: x(k+1)

1 se obtiene en función de x(k)2 , x

(k)3 , . . . , x

(k)n ; x(k+1)

2 se obtiene en funciónde x(k+1)

1 , x(k)3 , . . . , x

(k)n ; x(k+1)

3 en función de x(k+1)1 , x

(k+1)2 , x

(k)4 , . . . , x

(k)n ;. . . ; x(k+1)

n−1 en funciónde x

(k+1)1 , x

(k+1)2 , . . . , x

(k+1)n−2 , x

(k)n ; y x

(k+1)n en función de x

(k+1)1 , x

(k+1)2 , . . . , x

(k+1)n−1 . El hecho de

que el método de Gauss-Seidel utilice las componentes de la nueva aproximación �quenormalmente son más precisas que las de la aproximación anterior� antes que el método deJacobi, hace que �cuando converjan ambos métodos� el método de Gauss-Seidel converjacon más rapidez. Aun así, conviene conocer también el método de Jacobi para los casos enque este método converja y el método de Gauss-Seidel no converja.

El siguiente teorema determina cuándo converge el método de Gauss-Seidel.

Teorema 3.15 Sea A una matriz real inversible de orden n tal que aii = 0 para todo i =

1, 2, . . . , n; y sea b una matriz columna real del mismo orden. Consideramos la descomposi-

ción de A dada por (3.8) y la matriz GA de�nida por (3.13). Y para cualquier x(0) ∈ Rn, sea

{x(k)} la sucesión de puntos de Rn obtenida mediante la iteración (3.12) o (3.16). Entonces{x(k)} converge a la solución del sistema Ax = b si y solo si ρ(GA) < 1.

Los dos siguientes resultados proporcionan condiciones su�cientes (aunque no necesarias,a diferencia del teorema 3.15) para que converja el método de Gauss-Seidel. Son muy útilesen bastantes ocasiones. El primero de ellos tiene las mismas hipótesis que el teorema 3.13.

Teorema 3.16 Sea A una matriz real de orden n, y sea b una matriz columna real del

mismo orden. Si la matriz A es diagonalmente dominante (por �las o por columnas) en

sentido estricto, entonces la iteración de Gauss-Seidel (3.12) o (3.16) converge a la única

solución del sistema Ax = b para cualquier punto de partida x(0) ∈ Rn.

El segundo resultado exige una hipótesis menos que las que requería el teorema 3.13)para la convergencia del método de Jacobi.

62

Teorema 3.17 Si A es una matriz simétrica real de�nida positiva de orden n, y b es una

matriz columna real del mismo orden, entonces la iteración de Gauss-Seidel (3.12) o (3.16)converge a la única solución del sistema Ax = b para cualquier punto de partida x(0) ∈ Rn.

En muchas ocasiones ocurre que los métodos de Jacobi y de Gauss-Seidel no convergenpara un sistema dado pero, sin embargo, sí convergen para un sistema equivalente (obtenido,por ejemplo, permutando las ecuaciones o las incógnitas del sistema). Ilustraremos este hechoen los ejercicios.

3.6. Método de las potencias

En muchas aplicaciones del Álgebra Lineal es necesario encontrar el valor propio estric-tamente dominante (supuesto que exista) de una matriz cuadrada real. El radio espectralde la matriz coincide entonces con el valor absoluto10 de dicho valor propio. Por ejemplo, enla sección anterior nos hemos encontrado con problemas en los que fue necesario obtener elradio espectral de una matriz o, al menos, asegurar que dicho radio era menor que 1.

Uno de los métodos numéricos que proporcionan aproximaciones del valor propio estric-tamente dominante de una matriz cuadrada real es el llamado método de las potencias, quese expone a continuación.

Para que a una matriz real A de orden n (que se puede suponer no nula: sabemos quelas matrices nulas tienen como único valor propio λ = 0, con multiplicidad n) se le puedaaplicar el método de las potencias, debe cumplir las siguientes condiciones:

1. A tiene un valor propio estrictamente dominante, o bien un único valor propio domi-nante con multiplicidad mayor que 1 (que denotaremos por λ1);

2. A es diagonalizable.

Estas condiciones no son muy exigentes: la mayor parte de las matrices reales las cumplen.Tienen el inconveniente de que muchas veces no se pueden comprobar. En las prácticas conMathematica veremos cómo sortear este problema.

Supuesto que se cumplen ambas condiciones, partiendo de x(0) ∈ Rn − {θ}, se construyela siguiente sucesión de vectores de Rn:

x(k+1) = Ax(k), k = 0, 1, 2, . . .

Si existe el límite

lımk→∞

x(k+1)i

x(k)i

(3.17)

para algún i = 1, 2, . . . , n, entonces, bajo ciertas condiciones, ese límite será el valor propioλ1 buscado. Dichas condiciones son las siguientes:

1. Sea B = {v(1), v(2), . . . , v(n)} una base de Rn formada por vectores propios de la matrizA,11 donde v(1) es un vector propio no nulo asociado a λ1. Sean αi, i = 1, 2, . . . , n, lascoordenadas de x(0) con respecto a la base B (x(0) = α1v

(1) +α2v(2) + · · ·αnv

(n)). Puesbien, la condición es que se debe elegir x(0) de manera que α1 = 0.

10En la sección 3.1 vimos que los valores propios estrictamente dominantes son números reales.11El teorema 3.4 asegura la existencia de esa base.

63

2. El límite (3.17) debe tomarse para cualquier subíndice i tal que la componente i-ésimadel vector v(1) sea no nula: v(1)i = 0.

En general, estas condiciones no se podrán comprobar directamente. Para asegurar quese cumplen se recomienda calcular todos los límites (3.17) (para todo i = 1, 2, . . . , n), ypartir de distintos vectores iniciales x(0) que sean linealmente independientes entre sí. Así loharemos en las prácticas con Mathematica.

3.7. Ejercicios

1. Represente los siguientes conjuntos:

a) C2(ε) = {u ∈ R2 : ∥u∥2 ≤ ε};b) C1(ε) = {u ∈ R2 : ∥u∥1 ≤ ε};c) C∞(ε) = {u ∈ R2 : ∥u∥∞ ≤ ε}.

2. Obtenga los valores propios y el radio espectral de las siguientes matrices.

a) A =

(1 1

212

13

), E =

3 −1 1

0 2 1

0 1 2

, F =

3 0 0

0 2 0

0 1 2

, G =

3 0 −2

0 3 4

0 0 −3

.

b) [Mathematica] B =

1 0 1

−1 1 1

0 1 0

, C =

1 0 112

4 12

0 1 −32

, D =

1 0 112

4 12

−2 1 −32

.

3. a) Obtenga las normas matriciales ∥ · ∥∞ y ∥ · ∥1 de las matrices del ejercicio 2.

b) Obtenga también la norma ∥ · ∥2 de las matrices A, B, E, F y G de ese ejercicio.

c) Utilizando los dos apartados anteriores, acote el radio espectral de las siete ma-trices consideradas, y compare la cota con el valor del radio espectral (que se haobtenido en el ejercicio 2).

4. [Mathematica] Obtenga el condicionamiento (con respecto a la norma matricial quedesee) de la matriz de coe�cientes del siguiente sistema lineal:

−x1 + 2x2 + x3 = 2,

150x1 − 300x2 − 151x3 = −301,

250x1 − 500,1x2 − 250,9x3 = −501.

Resuelva de manera exacta el sistema anterior. Resuelva de nuevo el sistema si el vectorde los términos independientes fuera (2,01,−301,01,−500,99). Compare las solucionesde ambos sistemas. Obtenga el error relativo rx que se produce en la solución delsistema y compárelo con el error relativo en el vector de los términos independientesrb. Compare este error con las cotas dadas por la fórmula:

rbcondA

≤ rx ≤ condA · rb.

64

5. Repita el ejercicio anterior reemplazando la matriz de coe�cientes del sistema por lasiguiente: −1 2 1

0 3 1

−1 0 2

.

6. Obtenga el condicionamiento (con respecto a la norma matricial del in�nito) de lamatriz

1 0 2 −1

−1 0 0 2

1 4 2 0

0 0 −3 −1

.

7. Utilice los métodos de Jacobi y de Gauss-Seidel para resolver aproximadamente elsiguiente sistema de ecuaciones:

9x1 + x2 + x3 = 16,

2x1 + 10x2 + 3x3 = 44,

3x1 + 4x2 + 11x3 = 59.

a) Realice a mano tres iteraciones de cada método, comenzando desde el punto

(0, 0, 0). Acote el error cometido en cada caso.

b) [Mathematica] Utilice Mathematica hasta completar un total de 50 iteraciones.

8. Dado el sistema 3 1 7

5 1 0

1 4 2

x1

x2

x3

=

75

54

43

,

compruebe si existe alguna permutación de las �las que permita asegurar que los méto-dos de Jacobi y de Gauss-Seidel convergen.

9. Se considera el siguiente sistema de ecuaciones lineales:1 3 −1 0

1 2 3 −7

5 1 2 0

0 −1 3 1

x1

x2

x3

x4

=

1

−2

4

−1

.

a) Compruebe que existe alguna permutación de las ecuaciones del sistema que per-mite asegurar la convergencia del método de Jacobi.

b) Realice dos iteraciones de ese método, partiendo del punto (1, 1, 1, 1).

10. Se considera el siguiente sistema de ecuaciones lineales:2 6 −2 1

0 −1 3 1

2 3 4 −10

4 1 2 0

x1

x2

x3

x4

=

1

1

0

2

.

65

a) Compruebe que existe alguna permutación de las ecuaciones del sistema que per-mite asegurar la convergencia del método de Gauss-Seidel.

b) Pruebe que la fórmula que de�ne dicho método es

x(k) =

0 −1

4−1

20

0 112

12

−16

0 136

16

− 718

0 − 172

760

− 37180

x(k−1) +

12

0

13

730

, k ∈ N.

c) Partiendo del punto x(0) = (1/2, 0, 1/2, 0), aplique la iteración de Gauss-Seideluna vez; y encuentre una cota de la norma del error x(1) − x (use la norma quepre�era entre las que se pueden utilizar).

11. a) Dado el sistema 10 1 2 3

1 9 −1 2

2 −1 7 3

3 2 3 12

x1

x2

x3

x4

=

12

−27

14

−17

,

estudie si la matriz de coe�cientes es diagonalmente dominante (por �las o colum-nas) en sentido estricto y si es de�nida positiva.

b) Compruebe que su solución exacta es

(x1, x2, x3, x4) =

(9037

5559,−4704

1853,11038

5559,−3514

1853

).

c) Encuentre una aproximación de la solución partiendo de un punto próximo a esta(por ejemplo, x(0) = (1,5,−2,5, 2,−2)), realizando dos iteraciones de cada uno delos métodos que conozca cuya convergencia esté asegurada.

d) [Mathematica] Observe cuál de los dos métodos converge con mayor rapidez a lasolución.

12. Dadas las matrices

A =

3 2 2 −1

2 6 −1 2

2 −1 8 3

−1 2 3 7

y B =

5 −3 −3 0

−3 2 −2 1

−3 −2 4 2

0 1 2 7

,

pruebe que (a) los métodos de Jacobi y de Gauss-Seidel convergen para todo sistemalineal cuya matriz de coe�cientes es la matriz A, y que (b) [Mathematica] esos métodosno convergen para todo sistema lineal cuya matriz de coe�cientes es la matriz B.

66

13. Dado el sistema a b 0

1 3 1

0 a a

x1

x2

x3

=

a− b

0

b− a

,

determine una relación entre los números reales a y b para que el método de Jacobiconverja, y otra para que sea el método de Gauss-Seidel el que converja. Indicación:primero debe imponer que el sistema sea compatible determinado.

14. Determine una permutación de ecuaciones y una relación entre los valores de a y b, demodo que para el sistema a b a

a a b

b a a

x1

x2

x3

=

a− b

0

b− a

se tenga que: (a) el método de Jacobi converja; (b) el método de Gauss-Seidel converja.Estudie también cuándo la matriz del sistema permutado A = D − L− U es de�nidapositiva; y estudie también cuándo lo es la matriz D + L + U . ¾Qué conclusiones sepueden sacar de esto último?

15. [Mathematica] Consideramos el sistema 3 1 −1

0 3 1

−2 1 −4

x

y

z

=

1

3

2,5

.

a) Obtenga sendas aproximaciones de la solución del sistema, efectuando el númerode iteraciones de los métodos de Jacobi y de Gauss-Seidel que estime su�ciente.

b) Aplique de nuevo esos métodos si el coe�ciente que multiplica a la incógnita x3

en la tercera ecuación fuese −1 en lugar de −4.

c) Comente los resultados de los dos apartados anteriores, concretamente en cuantoa la convergencia de los métodos.

16. [Mathematica] Partiendo de la aproximación inicial x(0) = (0, 0, 0) (y de cualquier otraque elija), efectúe 40 iteraciones (extiéndalas después hasta 100 y 200) de los métodosde Jacobi y de Gauss-Seidel para aproximar �si fuera posible� la solución del sistema 15 17 19

0,25 0,33 0,42

1 1,2 1,6

x

y

z

=

2120

43,4

164

.

Muestre, en cada iteración, una cota del error de aproximación.

17. Realice 4 iteraciones del método de las potencias para aproximar el valor propio estric-tamente dominante de las siguientes matrices:

67

A =

(1 1

2

12

13

); B =

1 5 −8

5 −2 5

−8 5 1

; C =

1 2 3

2 3 4

3 4 5

;

D =

33 16 72

−24 −10 −57

−8 −4 −17

; E =

1 1

213

12

13

14

13

14

15

; F =

6 4 4 1

4 6 1 4

4 1 6 4

1 4 4 6

.

18. [Mathematica] La matriz de coe�cientes de cierto sistema es

2 −1 0

−1 2 −1

0 −1 2

.

Si es posible, utilice el método de las potencias para decidir si los métodos de Jacobi yde Gauss-Seidel son convergentes. Indicación: se trata de que el método de las potenciaspermita obtener aproximadamente el radio espectral de la matriz de cada método.

19. [Mathematica] Aplique el método de las potencias para obtener una aproximación,con al menos seis dígitos de precisión, del valor propio estrictamente dominante de lamatriz

A =

3 −3,5 −5,8 0 2

1 0 −3 2 2

−2 3 5 −0,5 0

−5 0 5 1 1

2 0 5 5 0

.

20. [Mathematica] Aplique el método de las potencias y realice el número de iteracionesnecesario para obtener una aproximación bastante precisa del valor propio estricta-mente dominante de la matriz:

A =

3,9375 −3,1875 −5,8125 0,25

1,875 −3,375 −5,625 2,5

−2,3125 3,0625 8,9375 −0,75

−5,625 5,125 16,875 1,5

.

21. [Mathematica] Aplique el método de las potencias para comprobar si el método deJacobi, aplicado a cualquier sistema lineal con matriz de coe�cientes 5 2 −3

2 7 2

0 2 1

es convergente o no. Repita el ejercicio para el método de Gauss-Seidel.

68

22. [Mathematica] Aplique el método de las potencias para comprobar si el método deGauss-Seidel, aplicado a cualquier sistema lineal con matriz de coe�cientes

5 0 −1 3

−1 5 2 −1

0 2 4 −1

3 −2 −1 −7

es convergente o no. Repita el ejercicio para el método de Jacobi. ¾Se puede concluiralgo acerca de la convergencia del método de Jacobi?

23. Se considera el siguiente sistema de ecuaciones lineales: 0 4 7

3 1 3

2 6 3

x

y

z

=

1

4

5

.

a) Compruebe que existe alguna permutación de las ecuaciones del sistema que per-mite asegurar la convergencia de los métodos de Jacobi y de Gauss-Seidel.

b) Una vez permutado adecuadamente el sistema, obtenga tres iteraciones del méto-do de Jacobi partiendo del punto que le parezca mejor.

24. Se considera el siguiente sistema de ecuaciones lineales:−1/2 −5 2 1

1/5 1 1 −7

−1 2 −3 4

1/4 1 −7 1

x

y

z

w

=

1

1

1

−1

.

Compruebe que existe alguna permutación de las ecuaciones del sistema que permiteasegurar la convergencia de los métodos de Jacobi y de Gauss-Seidel.

25. Dadas las matrices

A =

3 −1 0 1

0 3 −1 1

−1 3 2 −7

5 −6 −1 −2

y B =

(3 −1

5 2

),

obtenga ∥B∥2, ∥A∥1 y ∥A∥∞.

26. Obtenga una cota superior del error relativo de la solución del sistema

x+ 4y − z = 2,02 ,

2x− 3z = −3,005 ,

3x− 5y = −4,99 ,

si sus términos independientes se redondean a los valores 2, −3 y −5, respectivamente.

69

27. Se sabe que el método de Gauss-Seidel converge para el siguiente sistema:

4x− y + z = 1,

x− 4y + 3z = 0,

−3x− y + 5z = 3.

Partiendo de cierto punto, la iteración de Gauss-Seidel produce las siguientes aproxi-maciones de la solución del sistema:

(x4, y4, z4) = (0,1937, 0,6851, 0,8533) y (x5, y5, z5) = (0,2080, 0,6920, 0,8632).

Obtenga una cota de la norma del error cometido por esta última. Use la norma quepre�era.

28. Dado el sistema −1 9 −4

1 −1 2

2 −4 9

x

y

z

=

7

4

3

,

compruebe que existe una permutación de sus ecuaciones de manera que el método deGauss-Seidel �aplicado al sistema obtenido mediante dicha permutación� converja ala solución del sistema.

29. Dada la matriz

A =

−3 2 0

0 10 2

2 1 14

,

se desea aplicar el método de las potencias, partiendo del vector x(0) = (1, 1, 1), paracalcular aproximadamente el valor propio estrictamente dominante λ1 de A. Realicetres iteraciones de dicho método e indique las aproximaciones de λ1 que proporciona.

Algunas cuestiones teóricas

1. Sea A una matriz cuadrada real.

a) ¾Qué es el radio espectral de A?

b) ¾Qué relación hay entre el radio espectral de A y cualquiera de sus normas ma-triciales?

2. Sea A una matriz cuadrada real.

a) ¾Tiene A que tener obligatoriamente un valor propio estrictamente dominante?¾Si lo tiene, puede ser complejo? ¾Por qué?

b) Si A tiene un valor propio dominante real, ¾debe ser estrictamente dominante?¾Puede tener otro valor propio dominante que sea complejo? Razone las respues-tas.

70

3. En la teoría hemos estudiado tres distintas condiciones su�cientes para que el métodode Gauss-Seidel converja. Tomando dos de esas condiciones su�cientes, encuentre dossistemas compatibles determinados sencillos (de orden 2 o 3) en los que se cumpla unade las condiciones su�cientes pero no la otra, de manera que las condiciones que secumplen en cada sistema sean diferentes.

4. Dado un sistema compatible determinado Ax = b, donde A es una matriz simétricapero no de�nida positiva, ¾puede ocurrir que el método de Gauss-Seidel converja?Razone la respuesta.

5. Sea Ax = b un sistema lineal de orden n compatible determinado. Sea C la matriz deJacobi, la de Gauss-Seidel o la de cualquier otro método iterativo que converja a lasolución de dicho sistema (esto es, ρ(C) < 1). Sea x(k+1) = Cx(k) + d (k = 0, 1, 2, . . .)la iteración que de�ne el método, partiendo de cualquier x(0) ∈ Rn. Sea x la solucióndel sistema.

a) ¾Qué dos condiciones deben satisfacer las normas ∥ · ∥M y ∥ · ∥V para que puedaaplicarse la fórmula

∥x(k) − x∥V ≤ ∥C∥M1− ∥C∥M

∥x(k) − x(k−1)∥V , k = 1, 2, . . .?

b) En tal caso, explique qué utilidad tiene esa fórmula.

6. Indique al menos dos situaciones en las que el método de las potencias no converge.

7. Explique qué problema puede surgir al resolver un sistema lineal n × n compatibledeterminado si el condicionamiento de la matriz de coe�cientes es un número grande(por ejemplo, del orden de 106).

8. Supongamos que, para cierta matriz A se tiene que ∥A∥1 = 0,98, ∥A∥2 = 1,03 y∥A∥∞ = 0,92. ¾Qué se sabe entonces acerca del valor o valores propios dominantes deA? Razone la respuesta.

9. El radio espectral de una matriz A es igual a 4. (a) ¾Qué se sabe entonces acerca delos valores propios de A? (b) ¾Y de las normas matriciales de A?

10. ¾Es posible que el método de Gauss-Seidel converja para una matriz de la forma

A =

(a b

b 0

),

con a y b números reales no nulos? Razone la respuesta.

71

Capítulo 4

INTERPOLACIÓN DE FUNCIONES

La interpolación es un modo de aproximar funciones a partir de ciertos valores conocidosde ellas (y quizá también de sus derivadas). En este capítulo solo estudiaremos dos de losmétodos de interpolación más usados para funciones de una variable.1

4.1. Problemas de interpolación y de aproximación de

funciones

La interpolación de unos puntos dados de una función f consiste en encontrar otra funcióng, con ciertas características establecidas, que coincida con f en todos esos puntos. Es decir,si los puntos conocidos son los (xi, f(xi)), con i = 1, 2, . . . , n, se trata de encontrar unafunción g del tipo establecido tal que g(xi) = f(xi) para todo i = 1, 2, . . . , n. Entonces sedice que g es la función interpoladora y que g interpola a f en los puntos x1, x2, . . . , xn.

En ocasiones no se hace referencia a la función f , sino solo a unos datos o puntos ainterpolar (xi, yi), i = 1, 2, . . . , n, pero el problema es esencialmente el mismo; en este casose dice que g interpola los puntos (xi, yi) (i = 1, 2, . . . , n).

Un ejemplo práctico de este problema lo encontramos en el ejercicio 8 de la sección 4.4, enel que se tienen siete datos experimentales de la función �distancia recorrida por un automóvilque marcha a cierta velocidad hasta que se para (tras pisar el freno a fondo)�. Es decir, lafunción de partida f asigna a cada valor de la velocidad v la distancia recorrida tras frenarf(v), y esta función solo se conoce en siete puntos. El ejercicio pide que se obtenga unafunción polinómica g(v), de grado menor o igual que 6, que interpole esos siete puntos, esdecir, que su grá�ca pase por ellos. Cuando se obtenga g, se conocerá su valor para todovalor de v (a diferencia de la función de partida f).

También realizaremos ejercicios en los que la función f a interpolar se conozca en todosu dominio y, aun así, interese encontrar una función de otro tipo que la interpole en ciertospuntos. Nos encontraremos abundantes situaciones de este tipo en el capítulo 5. Pero el casomás frecuente en la práctica es el considerado en el párrafo anterior, en el que solo se conocenun conjunto �nito de datos de una función.

1La interpolación de funciones de dos o más variables, que tiene numerosas aplicaciones, escapa de las

posibilidades de este curso. Una introducción a la interpolación en dos variables puede encontrase en [11].

73

En algunos problemas de interpolación se dispone un mayor número de datos sobre cadapunto: por ejemplo, puede que se conozca no solo el valor de la función, sino también el desu primera derivada (e incluso más derivadas). En este caso se suele exigir que la funcióninterpoladora también coincida con la función de partida en esos valores de la derivada. Porejemplo, si se conocen los datos (xi, f(xi), f

′(xi)), con i = 1, 2, . . . , n, puede pedirse que lafunción interpoladora g cumpla que g(xi) = f(xi) y g′(xi) = f ′(xi) para todo i = 1, 2, . . . , n.En esta asignatura no estudiaremos este tipo de problemas de interpolación.

En cualquier caso, para cada problema práctico es importante elegir bien el tipo defunción interpoladora: su grá�ca debe asimilarse a la de la función a interpolar. Por ejemplo,cuando resolvamos el ejercicio 8 ya mencionado en clase de prácticas, se comprobará que esun error tomar un polinomio de grado menor o igual que 6 como función interpoladora g(v),puesto que produce una función que no es siempre creciente (y es evidente que la funcióninterpolada f(v) tiene que serlo: a mayor velocidad, más distancia recorrerá el automóvilhasta detenerse).

Hay otros muchos problemas en los que se busca aproximar una función f por otra funcióng de cierto conjunto G, de manera que g sea la función de G que se aproxima mejor �encierto sentido� a la función f . La interpolación puede entenderse como un problema deaproximación donde se exige que la función interpoladora se aproxime lo más posible a lafunción dada en los puntos a interpolar (en este caso se aproxima del todo, puesto que ambasfunciones son iguales en esos puntos). Pero hay muchos modos de entender que una funciónesté lo más próxima posible a otra. Dos de estos modos, quizá los más usados, son aquellosen los que la cercanía entre las funciones f y g se mide mediante las siguientes expresiones:

max{|f(x)− g(x)| : x ∈ D}, donde D es el dominio donde se de�nen f y las funcionesdel conjunto G;

n∑i=1

(f(xi)− g(xi))2, donde x1, x2, . . . , xn ∈ D y D es el mismo del caso anterior.

En ambos casos se busca la función g ∈ G que haga mínima la correspondiente expresión.Cuando la cercanía o distancia entre dos funciones se mide de la primera forma, tenemoslo que se llama un problema de aproximación uniforme; mientras que si se usa la segundaexpresión, es un problema de aproximación por mínimos cuadrados.

El lector conoce que, cuando se tiene un conjunto de datos {(xi, yi) : i = 1, 2, . . . , n} dedos variables estadísticas X e Y , la recta de regresión g0(x) de Y sobre X es precisamentela recta que hace mínima la expresión

n∑i=1

(yi − g(xi))2

entre todas las funciones lineales g(x) = ax+ b. En de�nitiva, los problemas de funciones deregresión son un caso particular de aproximación por mínimos cuadrados.

Los problemas de aproximación de funciones son tan interesantes o más que los proble-mas de interpolación, pero no hay espacio para ellos en esta asignatura. Una introduccióna la aproximación por mínimos cuadrados puede encontrarse en cualquier libro de teoríade la bibliografía; y en varios de ellos se puede encontrar también una introducción a laaproximación uniforme.

74

4.2. Interpolación polinomial: la fórmula de Newton

Los primeros métodos de interpolación utilizaban, por lo general, funciones polinómicascomo funciones interpoladoras. En esta sección introducimos uno de esos métodos.

Partimos de la situación descrita en la sección 4.1: se conocen los valores de una función f

en n puntos distintos y se desean interpolar estos datos. Buscaremos la función interpoladoraen el espacio vectorial Pn−1[x] de las funciones polinómicas de grado menor o igual que n−1.Para explicar cómo haremos esto, necesitamos introducir antes el siguiente concepto.

De�nición 4.1 Dada una función f y n puntos distintos de su dominio x1, x2, . . . , xn,2 sede�ne la diferencia dividida de la función f en esos puntos, y se denota f [x1, x2, . . . , xn], comoel número real de�nido mediante el siguiente proceso: f [xk] = f(xk) para todo k = 1, 2, . . . , n,y

f [xr, xr+1, . . . , xr+s] =f [xr, xr+1, . . . , xr+s−1]− f [xr+1, xr+2, . . . , xr+s]

xr − xr+s

(4.1)

para cualesquiera r, s ∈ N tales que 1 ≤ r < r + s ≤ n.

Así, para calcular la diferencia dividida f [x1, x2, . . . , xn], elaboraremos la siguiente tablade diferencias divididas:

x1 f [x1] f [x1, x2] f [x1, x2, x3] · · · f [x1, . . . , xn−2] f [x1, . . . , xn−1] f [x1, . . . , xn]

x2 f [x2] f [x2, x3] f [x2, x3, x4] · · · f [x2, . . . , xn−1] f [x2, . . . , xn]

x3 f [x3] f [x3, x4] f [x3, x4, x5] · · · f [x3, . . . , xn]...

......

...xn−2 f [xn−2] f [xn−2, xn−1] f [xn−2, xn−1, xn]

xn−1 f [xn−1] f [xn−1, xn]

xn f [xn]

Observe que la diferencia dividida f [x1, x2, . . . , xn] se sitúa en la última columna de la tabla;y para calcularla �recuerde (4.1)� será necesario obtener con anterioridad las que aparecenen la penúltima columna, que a su vez se obtendrán a partir de las diferencias divididas dela antepenúltima columna, y así sucesivamente. Luego habrá que comenzar por calcularlas diferencias divididas de la tercera columna en función de los datos de las dos primerascolumnas, que son los datos de partida de que dispondremos. Así, aplicando la de�nición(4.1) se obtiene

f [xi, xi+1] =f [xi]− f [xi+1]

xi − xi+1

, para todo i = 1, 2, . . . , n− 1.

Del mismo modo se obtienen las diferencias divididas de la cuarta columna en función de lasde la tercera:

f [xi, xi+1, xi+2] =f [xi, xi+1]− f [xi+1, xi+2]

xi − xi+2

, para todo i = 1, 2, . . . , n− 2.

2Los puntos x1, x2, . . . , xn que aparecen como argumentos de una diferencia dividida pueden estar orde-

nados de cualquier forma.

75

Y así sucesivamente, se obtienen las diferencias divididas de las siguientes columnas hastallegar a la última, cuyo único valor es

f [x1, x2, . . . , xn] =f [x1, x2, . . . , xn−1]− f [x1, x2, . . . , xn]

x1 − xn

.

Una propiedad interesante que satisfacen las diferencias divididas es que el orden en queaparecen los puntos en una diferencia dividida no altera su valor, es decir f [x1, x2, . . . , xn]

es igual a cualquier otra diferencia dividida de f en la que se hayan reordenado los puntosx1, x2, . . . , xn.

En el siguiente ejemplo se obtienen algunas diferencias divididas de una función.

Ejemplo 4.1 Dada la función f(x) = x2 + cos(πx) − 2 se desea calcular f [−2, 0, 1, 3]. Enprimer lugar se tiene que f [−2] = 3, f [0] = −1, f [1] = −2 y f [3] = 6. Por tanto, lasdiferencias divididas con dos argumentos que necesitamos calcular son

f [−2, 0] =f [−2]− f [0]

−2− 0= −2, f [0, 1] =

f [0]− f [1]

0− 1= −1, f [1, 3] =

f [1]− f [3]

1− 3= 4.

De aquí se obtienen las dos diferencias divididas con tres argumentos que necesitamos:

f [−2, 0, 1] =f [−2, 0]− f [0, 1]

−2− 1=

1

3, f [0, 1, 3] =

f [0, 1]− f [1, 3]

0− 3=

5

3.

Finalmente,

f [−2, 0, 1, 3] =f [−2, 0, 1]− f [0, 1, 3]

−2− 3=

4

15.

Las diferencias divididas obtenidas se representan en la siguiente tabla.

−2 f [−2] = 3 f [−2, 0] = −2 f [−2, 0, 1] = 13

f [−2, 0, 1, 3] = 415

0 f [0] = −1 f [0, 1] = −1 f [0, 1, 3] = 53

1 f [1] = −2 f [1, 3] = 4

3 f [3] = 6

Como se ha dicho, cualquier permutación de los argumentos da lugar al mismo resultado.Por ejemplo, f [0, 1,−2, 3] = f [−2, 0, 1, 3]. 2

El siguiente resultado asegura que el tipo de interpolación que vamos a estudiar en estasección produce una única función interpoladora del tipo elegido.

Teorema 4.1 Sea f una función real de�nida en un intervalo I. Supongamos que se conocen

los valores f(xi), con i = 1, 2, . . . , n, para n puntos distintos xi ∈ I. Entonces existe un

único polinomio p ∈ Pn−1[x] que interpola a f en esos puntos, es decir, que veri�ca que

p(xi) = f(xi) para todo i = 1, 2, . . . , n.

Aunque no demostramos el teorema, conviene observar que se veri�ca trivialmente paralos casos n = 1 y n = 2, ya que por un punto del plano solo pasa una recta horizontal, y pordos puntos del plano con distinta abscisa pasa una única recta que puede ser de cualquiertipo excepto vertical. En concreto, cuando se conoce un punto (x1, f(x1)) de la grá�ca def , entonces la única función polinómica constante p ∈ P0[x] que interpola a f en x1 es

76

p(x) = f(x1) para todo x ∈ I; y cuando se conocen dos puntos (x1, f(x1)) y (x2, f(x2)) de lagrá�ca de f , entonces la única función polinómica p ∈ P1[x] que interpola a f en x1 y x2 es

p(x) =(f(x2)− f(x1))x+ x2f(x1)− x1f(x2)

x2 − x1

, x ∈ I.

Para n = 3 el teorema 4.1 nos informa �entre otras cosas� de que para cualesquiera trespuntos no alineados del plano con abscisas distintas, existe una única función polinómica deltipo p(x) = ax2 + bx+ c, con a = 0, cuya grá�ca pasa por esos tres puntos.

El polinomio interpolador cuya existencia y unicidad asegura el teorema 4.1 se puedeobtener mediante la fórmula que proporciona el siguiente resultado.

Teorema 4.2 Sea f una función real de�nida en un intervalo I. Supongamos que se conoce

el valor de f en n puntos distintos x1, x2, . . . , xn ∈ I. Entonces el único polinomio p ∈ Pn−1[x]

que interpola a f en esos puntos es el dado por la siguiente fórmula, denominada fórmulade Newton del polinomio de interpolación:

p(x) =n∑

i=1

f [x1, x2, . . . , xi]i−1∏j=1

(x− xj). (4.2)

En la fórmula (4.2), cuando i = 1 se entiende que0∏

j=1

(x− xj) = 1. Así, si desarrollamos

el segundo miembro de (4.2) se obtiene

p(x) = f [x1] + f [x1, x2](x− x1) + f [x1, x2, x3](x− x1)(x− x2) + · · ·++ f [x1, x2, . . . , xn](x− x1)(x− x2) · · · (x− xn−1).

Veamos un ejemplo de aplicación de la fórmula de Newton.

Ejemplo 4.2 Se desea calcular el polinomio p(x) de grado menor o igual que 3 que interpolaa la función f(x) = x2 + cos(πx) − 2 en los puntos −2, 0, 1, 3. En el ejemplo 4.1 se obtuvoque f [−2] = 3, f [−2, 0] = −2, f [−2, 0, 1] = 1

3y f [−2, 0, 1, 3] = 4

15. Luego

p(x) = 3− 2(x+ 2) +1

3(x+ 2)x+

4

15(x+ 2)x(x− 1) = −1− 28x

15+

3x2

5+

4x3

15. 2

Si, después de obtener el polinomio de interpolación p(x) mediante (4.2), se necesitaseobtener el único polinomio q ∈ Pn[x] que interpola a f en los puntos dados y en un nuevopunto xn+1 (conociéndose el valor f(xn+1)), entonces la fórmula de Newton asegura que

q(x) = p(x) + f [x1, x2, . . . , xn, xn+1]n∏

j=1

(x− xj).

Luego no sería necesario repetir todos los cálculos, sino que bastaría calcular la diferenciadividida f [x1, x2, . . . , xn+1]; para lo cual se partiría de la tabla de diferencias divididas cons-truida para calcular f [x1, x2, . . . , xn], se le añadiría al �nal de las dos primeras columnaslos datos xn+1 y f(xn+1) y se añadiría una nueva diagonal al �nal de la tabla, haciendo loscorrespondientes cálculos. Esto es lo que se hace en el siguiente ejemplo.

77

Ejemplo 4.3 Se desea calcular el polinomio q(x) de grado menor o igual que 4 que interpolaa la función f(x) = x2 + cos(πx)− 2 en los puntos −2,−1, 0, 1, 3. Como en el ejemplo 4.2 seobtuvo que el polinomio p(x) de grado menor o igual que 3 que interpola a f en los puntos−2, 0, 1, 3 es p(x) = −1− 28x

15+ 3x2

5+ 4x3

15, entonces

q(x) = −1− 28x

15+

3x2

5+

4x3

15+ f [−2, 0, 1, 3,−1](x+ 2)x(x− 1)(x− 3).

Para calcular f [−2, 0, 1, 3,−1] necesitamos completar la tabla construida en el ejemplo 4.1de la siguiente forma:

−2 f [−2] = 3 f [−2, 0] = −2 f [−2, 0, 1] = 13

f [−2, 0, 1, 3] = 415

f [−2, 0, 1, 3,−1] = 25

0 f [0] = −1 f [0, 1] = −1 f [0, 1, 3] = 53

f [0, 1, 3,−1] = 23

1 f [1] = −2 f [1, 3] = 4 f [1, 3,−1] = 1

3 f [3] = 6 f [3,−1] = 2

−1 f [−1] = −2

Luego q(x) = −1− 28x

15+

3x2

5+

4x3

15+

2

5(x+ 2)x(x− 1)(x− 3). Operando se obtiene

q(x) = −1 +8x

15− 7x2

5− 8x3

15+

2x4

5. 2

4.3. Interpolación mediante funciones spline

En esta sección desarrollamos otro modo de interpolar puntos o funciones mediante ciertotipo de funciones polinómicas a trozos, las funciones spline (que de�nimos más abajo).

Cuando los datos que se deben interpolar no son pocos, e incluso cuando lo son, elpolinomio de interpolación suele presentar excesivas oscilaciones; y esto, por lo general, nose ajusta bien a los problemas que aparecen en la realidad. Es el caso del polinomio deinterpolación del ejercicio 8 �al que ya nos hemos referido en la sección 4.1�, que no escreciente en el dominio considerado (y debería serlo) por las oscilaciones que presenta sugrá�ca. Este problema se podrá evitar en este caso (y en otros muchos) eligiendo comofunción interpoladora una del tipo que de�nimos a continuación.

De�nición 4.2 Sean r ≥ 1 y n ≥ 2 dos números naturales. Una función spline �o, sim-plemente, un spline� de grado r con nodos x1, x2, . . . , xn (ordenados de la siguiente forma:x1 < x2 < · · · < xn) es una función s de�nida en el intervalo [x1, xn] que satisface las doscondiciones siguientes:

1. La restricción de s a cada intervalo [xi, xi+1] (i = 1, 2, . . . , n − 1), que denotaremossi, es una función polinómica de grado menor o igual que r; y al menos una de esasrestricciones si es de grado exactamente igual a r.

2. s es una función de clase Cr−1 en [x1, xn].3

3Las funciones de clase C0 en un intervalo (caso r = 1) son, por de�nición, las funciones continuas en

dicho intervalo.

78

Así, por ejemplo, un spline de grado 1 con nodos x1 < x2 < · · · < xn es una función con-tinua s de�nida en el intervalo [x1, xn], cuya restricción a [xi, xi+1] es una función polinómicade grado menor o igual que 1 para todo i = 1, 2, . . . , n− 1. Es decir, su grá�ca está formadapor la unión de los segmentos Li (con i = 1, 2, . . . , n − 1), donde Li es el segmento cuyospuntos extremos son (xi, s(xi)) y (xi+1, s(xi+1)). Luego hay un único spline de grado 1 connodos x1 < x2 < · · · < xn y con ordenadas respectivas y1, y2, . . . , yn ∈ R. Esta función tieneel inconveniente de que no suele ser útil para modelizar fenómenos reales, en los que lasfunciones suelen ser derivables. Pero en algunas situaciones puede ser aplicable. De hecho,en el capítulo 5 obtendremos una aplicación de los splines de grado 1.

A diferencia del caso r = 1, cuando r ≥ 2 se puede probar que hay in�nitos splines degrado r que interpolan los puntos (xi, yi), con i = 1, 2, . . . , n y x1 < x2 < · · · < xn. Es más,el número de esos splines depende de r − 1 parámetros.

Descartados los splines de grado 1 por no ser derivables, los splines cuadráticos (es decir,los de grado 2) tampoco son ideales: al no poder exigirse que sean de clase C2, en los nodos sepueden producir cambios bruscos en la convexidad (o concavidad) de la función. Sin embargoesto ya no sucede en los de grado 3. Y es claro que cuanto menor sea el grado del spline,menor será el gasto computacional para obtenerlo. Estos motivos, entre otros, han hecho quelos splines cúbicos (esto es, los de grado 3) sean los más utilizados.

En interpolación, el modo más usual de emplear funciones spline de cualquier grado esel que se explica a continuación. Supongamos que se desea interpolar unos puntos (xi, yi),donde i = 1, 2, . . . , n y x1 < x2 < · · · < xn, mediante una función spline. Los puntos dadosmuchas veces serán puntos de la grá�ca de una función f , de manera que yi = f(xi) paratodo i = 1, 2, . . . , n.4 La situación más sencilla, que es la que abordaremos aquí, es cuando sedecide que los nodos del spline buscado coincidan con las abscisas de los puntos a interpolar.

En el siguiente ejemplo se encuentra un spline cuadrático que interpola unos puntos dadosde la manera que acabamos de explicar.

Ejemplo 4.4 Dados los puntos (−1, 2), (1, 0) y (2,−2), obtengamos los splines cuadráticosque los interpolan, de modo que los nodos de esos splines sean las abscisas de dichos puntos:los nodos son x1 = −1, x2 = 1 y x3 = 2. Por tanto, si s es uno de esos splines, s está de�nidoa trozos mediante dos funciones cuadráticas s1(x) = a1x

2+b1x+c1 y s2(x) = a2x2+b2x+c2,

de�nidas en los intervalos [−1, 1] y (1, 2], respectivamente.En primer lugar, para que s interpole los tres puntos dados y sea además una función

continua, tiene que suceder que s1(−1) = 2, s1(1) = 0 = lımx→1 s2(x) y s2(2) = −2. Es decir,deben satisfacerse las cuatro ecuaciones siguientes:

a1 − b1 + c1 = 2,

a1 + b1 + c1 = 0,

a2 + b2 + c2 = 0,

4a2 + 2b2 + c2 = −2.

Y para que s sea derivable debe suceder además que s′1(1) = lımx→1 s′2(x) (con esta

condición también se asegura ya que s es de clase C1 en el intervalo [−1, 2]). Como s′1(x) =

4En de�nitiva, seguimos en las condiciones generales comentadas en los párrafos 1 y 3 de la sección 4.1.

79

2a1x+ b1 y s′2(x) = 2a2x+ b2, se llega a la ecuación

2a1 + b1 = 2a2 + b2,

que completa con las cuatro anteriores un sistema lineal de 5 ecuaciones con 6 incógnitas,cuya solución es

(a1, b1, c1, a2, b2, c2) =(− c2

4,−1, 1 +

c24,c22− 1, 1− 3c2

2, c2

)para todo c2 ∈ R.

Es decir, hay in�nitos splines cuadráticos que interpolan los puntos dados, que dependen deun parámetro c2 ∈ R:

s(x) =

−c2

4x2 − x+ 1 +

c24

x ∈ [−1, 1],(c22− 1)x2 +

(1− 3c2

2

)x+ c2 x ∈ (1, 2],

Por ejemplo, para c2 = 0 se obtiene el spline

s(x) =

{1− x x ∈ [−1, 1],

x− x2 x ∈ (1, 2].

También podría determinarse un spline concreto imponiendo que pase por un puntodeterminado cuya abscisa no sea un nodo del spline (serían condiciones del tipo s(x0) = y0,con x0, y0 ∈ R y x0 = −1, 1, 2). 2

El siguiente teorema proporciona un método simple �y fácilmente programable� paraobtener el spline cúbico que interpola unos puntos dados de manera que las abscisas de estossean los nodos del spline.

Teorema 4.3 Supongamos que unos puntos dados (xi, yi), con i = 1, 2, . . . , n, veri�can que

x1 < x2 < · · · < xn. Entonces, una función s es un spline cúbico de nodos x1, x2, . . . , xn que

interpola a esos n puntos si, y solo si, las restricciones si de s a los intervalos [xi, xi+1], con

i = 1, 2, . . . , n− 1, son de la siguiente forma:

si(x) =1

6(xi+1 − xi)

[di+1(x− xi)

3 + di(xi+1 − x)3 +

+(6yi+1 − di+1(xi+1 − xi)

2)(x− xi) +

(6yi − di(xi+1 − xi)

2)(xi+1 − x)

],

donde (d1, d2, . . . , dn) es cualquier solución del sistema lineal formado por las n−2 ecuaciones

siguientes:

(xi − xi−1)di−1 + 2(xi+1 − xi−1)di + (xi+1 − xi)di+1 = 6

(yi+1 − yixi+1 − xi

− yi − yi−1

xi − xi−1

), (4.3)

con i = 2, 3, . . . , n − 1. Además, en tal caso resulta que cada di es la derivada segunda del

spline cúbico en el nodo xi: di = s′′(xi) para todo i = 1, 2, . . . , n.

80

A continuación estudiamos el sistema lineal a que se re�ere el teorema 4.3. Observe queen cada ecuación (4.3) aparecen exactamente tres incógnitas, que son consecutivas; y que laincógnita d1 solo aparece en la primera ecuación, mientras que la incógnita dn solo apareceen la última. Si en estas dos ecuaciones pasamos los sumandos que contienen a d1 o dn alsegundo miembro, entonces esas ecuaciones quedan así:

2(x3 − x1)d2 + (x3 − x2)d3 = 6

(y3 − y2x3 − x2

− y2 − y1x2 − x1

)− (x2 − x1)d1,

(xn−1 − xn−2)dn−2 + 2(xn − xn−2)dn−1 = 6

(yn − yn−1

xn − xn−1

− yn−1 − yn−2

xn−1 − xn−2

)− (xn − xn−1)dn,

donde d1 y dn pasan a considerarse parámetros reales en vez de incógnitas. Con esta premisa,el sistema ha pasado a ser un sistema de orden n − 2, y puede describirse matricialmentecomo Hd = b, donde

H =

2(x3 − x1) x3 − x2 0 · · · · · · 0

x3 − x2 2(x4 − x2) x4 − x3 0 · · · 0

0 x4 − x3. . . . . . . . .

...... 0

. . . . . . . . . 0...

.... . . . . . . . . xn−1 − xn−2

0 0 · · · 0 xn−1 − xn−2 2(xn − xn−2)

,

d =

d2d3...

dn−1

y b =

6

(y3 − y2x3 − x2

− y2 − y1x2 − x1

)− (x2 − x1)d1

6

(y4 − y3x4 − x3

− y3 − y2x3 − x2

)...

6

(yn−1 − yn−2

xn−1 − xn−2

− yn−2 − yn−3

xn−2 − xn−3

)6

(yn − yn−1

xn − xn−1

− yn−1 − yn−2

xn−1 − xn−2

)− (xn − xn−1)dn

.

Observe que el sistema Hd = b es tridiagonal (es decir, su matriz de coe�cientes escuadrada y, fuera de los elementos de la diagonal principal y de las dos diagonales situadasinmediatamente encima y debajo de la principal, todos los elementos son nulos). Puededemostrarse que el determinante de H es distinto de cero, por lo que d = H−1b expresa lasolución del sistema.

En de�nitiva, el sistema formado por las ecuaciones (4.3), con i = 2, 3, . . . , n − 1, tieneuna única solución para cada elección de los parámetros d1 y dn.

Cuando se establece que d1 = dn = 0, entonces se obtiene el denominado spline cúbico

natural.5

5Este spline cúbico se llama natural porque sería aquel cuya grá�ca podría prolongarse de forma �natural�

�es decir, de manera que se mantenga la suavidad de la curva� mediante una semirrecta por cada extremo.

En efecto, si se prolongara la de�nición de este spline fuera del intervalo [x1, xn] mediante dos funciones

lineales, de manera que existiera la derivada del spline �prolongado de esta manera� en los nodos x1 y

81

4.4. Ejercicios

1. Demuestre que la función g(x) = (eπ/2 + e−π/2) sen x + 2 cosx interpola a la funciónf(x) = ex + e−x en los puntos (0, 2) y (π/2, eπ/2 + e−π/2).

2. Demuestre que la función polinómica q(x) = 7π√3x − 3πx2 interpola a la función

f(x) = 36 arc tg(x) en los puntos (0, 0), (1/√3, 6π) y (

√3, 12π).

3. Pruebe que el polinomio p ∈ P4[x] que interpola los puntos (0, 5), (2, 13), (3, 2), (4,−1)

y (6, 10) es

p(x) = 5 +113x

3− 1963x2

72+

97x3

16− 61x4

144.

4. Obtenga el polinomio p ∈ P4[x] que interpola los puntos (−4, 0), (−1, 0), (1, 2), (3, 1)y (5, 3).

5. [Mathematica] Obtenga el polinomio p de grado no mayor que 7 que interpola a lafunción f(x) = e−x3/2 cos x en los puntos xi = i/10, con i = 1, 2, . . . , 8. Aplique lafórmula obtenida para dar una aproximación de f(0,26), y compárela con su verdaderovalor. Represente grá�camente los puntos dados y su polinomio interpolador.

6. [Mathematica] Repita el ejercicio 5 si los puntos donde se interpola la función sonx1 = 0, x2 = 2, x3 = 1,5, x4 = 1,3, x5 = 0,4, x6 = 1,300001 y x7 = 0,2.

7. Se sabe que la ecuación x − 9−x = 0 tiene una solución en el intervalo [0, 1]. Calculeel polinomio de interpolación de la función del primer miembro de esa ecuación en lospuntos 0, 0,5 y 1, y utilícelo para dar una solución aproximada de dicha ecuación.

8. [Mathematica] La distancia requerida para detener un automóvil que marcha a ciertavelocidad tras pisar el freno (a fondo) es una función de dicha velocidad. Se han recogidolos siguientes datos experimentales para cuanti�car esa relación:

velocidad (km/h) 35 40 50 55 60 70 80

distancia (m) 4,8 6 10,2 12 18 27 36

a) Utilice el polinomio que interpola estos datos para estimar la distancia de frenadodel automóvil si frena a fondo cuando marcha a una velocidad de 65 km/h.

b) Teniendo en cuenta los resultados del apartado anterior, ¾le parece acertado hacerestimaciones de la distancia de frenado usando el polinomio de interpolación? ¾Porqué?

9. Obtenga el polinomio p de grado menor o igual que 2, que tome los valores 1, 2 y −1

en los puntos 0, 1 y −2, respectivamente. Obtenga también el polinomio q ∈ P3[x] queinterpole los puntos anteriores y que cumpla que q(−1) = 1.

xn (con lo que la pendiente de la función lineal situada a la izquierda de x1 sería s′(x1), y la pendiente de

la función lineal situada a la derecha de xn sería s′(xn), entonces la función resultante seguiría siendo de

clase C2, esta vez en R, ya que la segunda derivada de una función lineal es cero, como lo son s′′(x1) = d1 y

s′′(xn) = dn, por hipótesis.

82

10. Halle los polinomios que interpolan a las siguientes funciones en los puntos que seindican.

a) f(x) =10

1 + x2en los puntos −2, 1 y 3.

b) f(x) = x3 en los puntos −1, 0, 1 y 3.

Dé una razón que justi�que el resultado que se obtiene en el apartado b).

11. a) Obtenga el polinomio que interpola a la función f(x) = x4 en los puntos −1, 1, 2.

b) Obtenga el polinomio que interpola a la función f(x) = x4 en los puntos −1, 1, 2y 3.

c) Compare los valores de los dos polinomios de interpolación obtenidos con los dela propia función f en los puntos −0,5, 0,5, 1,5 y 2,5. ¾Le parece aceptable lainterpolación realizada �en cuanto a la proximidad del polinomio a la funciónque interpola� en alguno de los dos casos?

12. Obtenga los polinomios que interpolan a las siguientes funciones en los puntos que seindican.

a) f(x) =1

2x+ 1, en los puntos −1, 0 y 1.

b) g(x) =1

1 + x2, en los puntos −2, −1, 1 y 2.

c) h(x) = cos2 x, en los puntos −0,5, 0 y 0,5.

13. Sea q(x) el polinomio que interpola los datos (0, 0), (1/2, b), (1, 3) y (2, 2), donde b ∈ R.¾Para qué valor de b se tiene que el grado de q es 2?

14. [Mathematica] Obtenga el polinomio p del grado más bajo posible que interpola lospuntos (0, 5), (1, 7), (2, 13), (3, 2), (4,−1), (5,−2) y (6, 13) de la grá�ca de ciertafunción f . ¾Cuál sería el valor aproximado de f(2,3) que proporciona el polinomio deinterpolación? Obtenga también el polinomio q que interpola los puntos anteriores ytambién el punto (7, 29). ¾Qué valor aproximado de f(2,3) proporciona q? Representelos dos polinomios en el mismo grá�co.

15. Obtenga el spline cúbico natural s(x) con nodos −1, 0, 1 tal que s(−1) = 13, s(0) = 7

y s(1) = 9.

16. [Mathematica] Obtenga la expresión general de los splines cúbicos cuyos nodos sonlos cinco primeros números naturales y que interpolan los puntos (1, 0), (2,−1), (3, 3),(4, 3) y (5, 1). Además, obtenga dos casos particulares:

a) El spline cúbico natural.

b) El spline cúbico s(x) tal que s′′(1) = −100 y s′′(5) = 100.

Finalmente, represente en un mismo grá�co los dos splines anteriores.

17. Obtenga los splines cúbicos que interpolan los puntos del ejercicio 4 y cuyos nodos sonprecisamente las abscisas de esos puntos.

83

18. a) Obtenga los splines cúbicos que interpolan a la función f(x) = x4 en los nodos−1, 1, 2 y 3. En particular, obtenga el spline cúbico natural y su valor en cada unode los puntos −0,5, 0,5, 1,5 y 2,5; y utilice los resultados obtenidos en el apartadoc) del ejercicio 11 para comparar esos valores del spline con los que toma f y elpolinomio que interpola a f en −1, 1, 2 y 3 en dichos puntos −0,5, 0,5, 1,5 y 2,5.

b) Obtenga el spline cúbico s que interpola a f en los nodos −1, 1, 2 y 3 y quetambién interpole a f en los puntos 0 y 1/3 (note que estos dos puntos no sonnodos del spline). Obtenga el valor del s en los puntos −0,5, 0,5, 1,5 y 2,5, ycompare estos valores con los considerados en el apartado a) de este ejercicio.

19. [Mathematica] Repita el ejercicio 8 interpolando con el spline cúbico natural. Comentesi el resultado le parece satisfactorio.

20. Compruebe si la siguiente función es un spline cuadrático con nodos 0, 1, 2 y 3.

h(x) =

x, x ∈ [0, 1],

3− (2− x)2

2, x ∈ (1, 2],

3/2, x ∈ (2, 3].

21. Determine para qué valores a, b, c, d ∈ R se cumple que la función

f(x) =

{3 + x− 9x2, x ∈ [0, 1],

a+ b(x− 1) + c(x− 1)2 + d(x− 1)3, x ∈ (1, 2]

es un spline cúbico con nodos 0, 1 y 2. ¾Puede ser f en algún caso un spline cúbiconatural?

22. Determine las condiciones que deben satisfacer los números reales a, b, c, d y e paraque la función

g(x) =

a+ bx− 39x2 + 13x3, x ∈ [1, 2],

232− 314x+ cx2 + dx3, x ∈ (2, 4],

−1624 + 1078x− 225x2 + ex3, x ∈ (4, 5]

sea un spline cúbico con nodos 1, 2, 4 y 5. Compruebe también si g puede ser un spline

cúbico natural.

23. [Mathematica] Dados los puntos (−4, 7), (−2, 5), (0, 13), (1, 7), (2, 1), (4,−1), (5,−2),obtenga los splines cúbicos y el polinomio de P6[x] que los interpole. Represente estepolinomio y el spline cúbico natural en el mismo grá�co.

24. Obtenga el polinomio que interpola los puntos (0, 0), (√3/3, π/6) y (

√3, π/3).

25. a) Obtenga el polinomio p que interpola a la función f(x) = ex en los puntos ln 2,ln 3 y ln 5.

b) ¾Qué habría que sumar al polinomio p para obtener el polinomio q que interpolaa f en los puntos anteriores y también en el punto 0?

84

26. a) Determine para qué valores reales de A, B, C, D y E la función

f(x) =

A(1 + x+ x2 + x3), x ∈ [a, 0],

B + x+ x2 + Cx3, x ∈ (0, 1],

D + 10x+ Ex2 + 2x3, x ∈ (1, b]

es un spline cúbico con nodos a, 0, 1 y b (donde a < 0 y b > 1).

b) ¾Es posible que el spline cúbico del apartado anterior sea natural para ciertosvalores de a y b?

27. Halle los números a, b, c, d ∈ R para que la función

f(x) =

{a(x− 1)3 + b(x− 1) + 2, 1 ≤ x < 2,

c(x− 2)3 − 3(x− 2)2 + d(x− 2) + 4, 2 ≤ x ≤ 3

sea un spline cúbico natural que interpola los puntos (1, 2), (2, 4) y (3, 2).

28. Dados los puntos (1, 2), (2, 0) y (4, 2), obtenga los splines cuadráticos que los interpolan,de modo que los nodos de esos splines sean las abscisas de dichos puntos. En particular,obtenga el spline cuadrático que satisface que su primer trozo es de grado 1.

29. Dados los puntos (−1, 0), (1, 2), (2, 6) y (3, 4), obtenga el spline cuadrático que losinterpola, de modo que los nodos de ese spline sean las abscisas de dichos puntos yel primer trozo del spline sea de grado 1. Indicación: Este ejercicio requiere muchasoperaciones.

Algunas cuestiones teóricas

1. Dados n, k ∈ N, ¾qué función polinómica interpola en n + k puntos distintos a otrafunción polinómica de grado n? Indicación: dé valores concretos a n y k y luego intentegeneralizar el resultado.

2. Compruebe que las funciones polinómicas q(x) = 72 − 6x − 55x2 − 15x3 + 3x4 + x5

y p(x) = 24 − 10x − 15x2 + x4 interpolan ambas los puntos (−3, 0), (−2, 0), (1, 0),(4, 0). ¾Contradice este hecho el teorema 4.1, donde se habla de un único polinomio deinterpolación?

3. Responda a las siguientes cuestiones razonando las respuestas y sin hacer cálculos.

a) Determine cuál es el polinomio de menor grado posible que interpola los siguientespuntos (1,−1), (−1, 1), (2,−2), (−2, 2), (3,−3), (−3, 3), (4,−4), (−4, 4), (5,−5)

y (−5, 5).

b) Teniendo en cuenta el apartado anterior, elija cuál de las siguientes frases es lacorrecta: (1) �Es POSIBLE (pero NO SEGURO) que exista un polinomio de grado9 que interpole dichos puntos�. (2) �Es IMPOSIBLE que exista un polinomio degrado 9 que interpole dichos puntos�. (3) Es SEGURO que existe un polinomiode grado 9 que interpola dichos puntos�.

85

Capítulo 5

INTEGRACIÓN NUMÉRICA

La regla de Barrow permite integrar de modo exacto las funciones cuya primitiva existey puede calcularse (aunque sea a trozos). Muchas veces esto no es posible, bien porque lafunción solo se conoce en una serie de puntos, o porque no hay modo de expresar la primitivade una forma elemental. En estos casos, que son muy frecuentes, se hace necesario encontrarmétodos numéricos que aproximen el valor de la integral deseada.

Los métodos de integración numérica que vamos a estudiar son de tipo interpolatorio.Esto quiere decir que la integral de una función f en un intervalo [a, b] se aproximará de lasiguiente forma: ∫ b

a

f(x) dx ≃∫ b

a

s(x) dx, (5.1)

donde s es una función que interpola a f en determinados puntos de [a, b] y que es fácil deintegrar en este intervalo.

5.1. Fórmulas simples

En esta primera sección establecemos las bases para obtener fórmulas de integraciónnumérica de tipo interpolatorio que proporcionen buenas aproximaciones de la integral exac-ta: esto no lo conseguiremos hasta la sección siguiente, pero basándonos en lo que estudiamosen esta. A las fórmulas que introducimos en esta sección se les denomina fórmulas simples deintegración numérica, por contraposición a las fórmulas compuestas de integración numéri-ca de la sección siguiente; estas últimas se obtendrán aplicando repetidamente las fórmulassimples.

A continuación estudiamos los métodos más sencillos de obtener una función s que in-terpole a una función dada f , ambas de�nidas en un intervalo [a, b], de manera que puedaestablecerse la aproximación (5.1). En estos métodos sencillos s será un polinomio de gradobajo que interpole a f en unos puntos según lo estudiado en la sección 4.2.

En primer lugar, consideraremos el caso en que s es el polinomio de P1[x] que interpolaa f en los extremos del intervalo, es decir,

s(x) = f(a) +f(b)− f(a)

b− a(x− a)

87

para todo x ∈ [a, b] (fórmula que puede obtenerse directamente o mediante la fórmula (4.2),con lo que la fórmula (5.1) adopta ahora la siguiente forma:∫ b

a

f(x) dx ≃ b− a

2(f(a) + f(b)). (5.2)

Esta es la llamada fórmula del trapecio.1

En segundo y último lugar, consideramos el caso en que s es el polinomio de P2[x] queinterpola a f en los extremos del intervalo y también en su punto medio. Aproximando denuevo la integral de f por la integral de s, se obtiene la siguiente fórmula:∫ b

a

f(x) dx ≃ b− a

6

[f(a) + 4f

(a+ b

2

)+ f(b)

]. (5.3)

Esta es la llamada fórmula de Simpson.

5.2. Fórmulas compuestas

En esta sección mostramos el modo de obtener fórmulas de integración numérica detipo interpolatorio que sean verdaderamente útiles para aproximar integrales de�nidas defunciones reales de variable real. El método no es otro que aplicar fórmulas simples en lossubintervalos en que se va a dividir el intervalo de integración. De ahí que a las fórmulasque vamos a introducir se denominan fórmulas compuestas de integración numérica (sonfórmulas �compuestas� de otras simples).

Supongamos que se desea aproximar el valor de la integral de una función integrable f

en un intervalo [a, b]. Se divide este intervalo en n partes iguales de longitud

h =b− a

n.

Por tanto, los puntos que determinan la partición de [a, b] son a = x0 < x1 < · · · < xn = b,con

xi = a+ ih para todo i = 0, 1, . . . , n.

Entonces, las fórmulas compuestas de integración numérica se obtendrán aplicando fórmulassimples a cada uno de los sumandos del segundo miembro de la siguiente expresión:∫ b

a

f(x) dx =n∑

i=1

∫ xi

xi−1

f(x) dx. (5.4)

En esta sección nos limitaremos a estudiar dos casos de fórmulas compuestas: las que seobtienen aplicando las fórmulas simples del trapecio y de Simpson, que proporcionan dos delos métodos de integración numérica más utilizados y también más sencillos.

1Se denomina así porque la región delimitada por s y el intervalo [a, b] en el eje de abscisas es un trapecio

cuando f(a) y f(b) son del mismo signo.

88

5.2.1. Método de los trapecios

El método de los trapecios consiste en aplicar la fórmula simple del trapecio en cadasumando de la expresión (5.4), de manera que se aproxima la integral de la función f en[a, b] por la integral de la función cuya grá�ca es la línea poligonal formada por cada uno delos segmentos que unen los puntos (xi−1, f(xi−1)) y (xi, f(xi)) (i = 1, 2, . . . , n): en la �gura5.1 se representan algunos de los trapecios asociados a una función y a la partición elegida.

O a x

y

b

f(x)

x

xx2

n-1

1...

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

f(b)f(x )

f(x )

f(x )f(a)

n-1

1

2

Figura 5.1: Representación del método de los trapecios

Es intuitivamente claro que la aproximación al valor exacto de la integral de f en [a, b]

tenderá a ser mejor a medida que el número n de puntos de la partición aumente. Esto locon�rmaremos en la sección 5.3.

Aplicando la fórmula (5.2) en cada sumando de (5.4) se obtiene∫ b

a

f(x) dx ≃n∑

i=1

xi − xi−1

2

(f(xi−1) + f(xi)

),

de donde se deduce inmediatamente la siguiente fórmula, conocida como fórmula compuesta

de los trapecios (o simplemente fórmula de los trapecios):∫ b

a

f(x) dx ≃ h

2

[f(a) + 2

n−1∑i=1

f(xi) + f(b)

]. (5.5)

5.2.2. Método de Simpson

El método de Simpson consiste en aplicar la fórmula simple de Simpson en cada sumandode la expresión (5.4), de manera que se aproxima la integral de la función f en [a, b] por laintegral de una función cuadrática a trozos (su grá�ca está formada por trozos de parábolas):en cada trozo [xi−1, xi] se integra el polinomio de grado menor o igual que 2 que interpola af en los puntos xi−1, xi− 1

2y xi, donde

xi− 12=

xi−1 + xi

2= a+

(i− 1

2

)· h para todo i = 1, 2, . . . , n. (5.6)

89

Aplicando la fórmula (5.3) en cada sumando de (5.4) se obtiene∫ b

a

f(x) dx ≃n∑

i=1

xi − xi−1

6

(f(xi−1) + 4f(xi− 1

2) + f(xi)

),

de donde se deduce inmediatamente la fórmula compuesta de Simpson (nos referiremos a ellasimplemente como la fórmula de Simpson, si no hay peligro de confundirla con la fórmulasimple del mismo nombre):∫ b

a

f(x) dx ≃ h

6

[f(a) + 2

n−1∑i=1

f(xi) + 4n∑

i=1

f(xi− 12) + f(b)

]. (5.7)

Observe que la fórmula (5.7) utiliza 2n+ 1 evaluaciones de la función f , mientras que lafórmula (5.5) solo utiliza n+1. Por tanto, el esfuerzo computacional de la fórmula de Simpsonque utiliza n subintervalos equivale al de la fórmula de los trapecios con 2n subintervalos.Visto al revés, que es lo más natural, si tenemos un número impar k de puntos equidistantesde la recta real donde se conoce el valor de una función f , para aproximar la integral de f

entre el menor y el mayor de dichos puntos se podrá utilizar la fórmula de los trapecios paran = k− 1 o la fórmula de Simpson para n = k−1

2. Si el número de puntos que se tiene es par,

solo puede utilizarse la fórmula de los trapecios.

5.3. Estudio del error de las fórmulas compuestas

Como en todo método numérico, es importante tener cierto conocimiento del error queproduce: es necesario saber hasta qué punto el valor aproximado dado por el método seacerca al valor exacto de la integral. A este punto dedicamos esta sección.

Intuitivamente, el modo normal de disminuir el error será aumentar n, es decir, el númerode subintervalos en que se divide el intervalo de integración o, lo que es lo mismo, aumentarel número de puntos en los que hay que evaluar la función del integrando para obtener lafórmula compuesta (n+1 en la fórmula de los trapecios y 2n+1 en la de Simpson). De hecho,el estudio que haremos del error cometido por los métodos de los trapecios y de Simpsonpermitirá encontrar el número de intervalos a los que aplicar la fórmula compuesta para quese asegure que el error de aproximación sea menor que una cantidad preestablecida.

5.3.1. Error de la fórmula de los trapecios

Bajo hipótesis bastante generales, el siguiente teorema proporciona una expresión delerror de la fórmula de los trapecios.

Teorema 5.1 Sea f una función de clase C2 en [a, b]. Para cierto n ∈ N, sea h = (b−a)/n

y sea a = x0 < x1 < · · · < xn = b la partición de [a, b] de�nida mediante xi = a + ih para

todo i = 0, 1, 2, . . . , n. Entonces existe un número c ∈ [a, b] tal que∫ b

a

f(x) dx =h

2

[f(a) + 2

n−1∑i=1

f(xi) + f(b)

]− (b− a)3

12 · n2f ′′(c).

90

Al ser f ′′ continua en [a, b], existe M > 0 tal que |f ′′(x)| ≤ M para todo x ∈ [a, b]. Loideal será tomar

M = max{|f ′′(x)| : x ∈ [a, b]}, 2

siempre que pueda determinarse dicho máximo. En de�nitiva, el error absoluto cometido conla fórmula de los trapecios se puede acotar, ya que∣∣∣∣− (b− a)3

12 · n2f ′′(c)

∣∣∣∣ = (b− a)3

12 · n2

∣∣f ′′(c)∣∣ ≤ M(b− a)3

12 · n2.

En el último término de la expresión anterior, M , b y a son cantidades �jas. Por tanto,para hacer la cota del error tan pequeña como se desee bastará aumentar el número n desubintervalos del método. En concreto, si se desea que el error absoluto sea menor que ciertacantidad positiva δ, bastará imponer que

M(b− a)3

12 · n2< δ,

es decir, n deberá satisfacer que

n >

√M(b− a)3

12 · δ=

b− a

2

√M(b− a)

3 · δ. (5.8)

Un ejemplo de aplicación de lo que acabamos de mostrar es el siguiente.

Ejemplo 5.1 Se desea aproximar la integral∫ 8/9

0

cosx2 dx

mediante el método de los trapecios con un error absoluto inferior a una milésima. Para locual se determina el número n de subintervalos que asegure que el error no superará dichacota. Luego n debe satisfacer (5.8) para δ = 10−3. Para determinar la constante M en (5.8)se calcula la segunda derivada de la función f(x) = cosx2, que es

f ′′(x) = −4x2 cosx2 − 2 sen x2.

Puede comprobarse que f ′′′(0) = 0 y que f ′′′(x) < 0 para todo x ∈ (0, 8/9], por lo que f ′′ esestrictamente decreciente en [0, 8/9]. Como f ′′(0) = 0 se tiene que

max{|f ′′(x)| : x ∈ [0, 8/9]} = |f ′′(8/9)| ≃ |−3,6451| = 3,6451.

Luego podemos tomar M = 3,6451 y, por tanto,

n >b− a

2

√M(b− a)

3 · δ≃ 4

9

√3,6451 · (8/9)

3 · 10−3=

8

27

√2 · 3645,1

3≃ 14,6061.

Es decir, si n ≥ 15 se puede asegurar que la integral se ha aproximado con un error inferiora 10−3. 2

2El teorema de los valores extremos asegura que este máximo existe.

91

En ocasiones se podrá obtener una aproximación algo mejor de la integral si se conocenuna cota inferior M1 y otra superior M2 de f ′′(x) en el intervalo [a, b]: M1 ≤ f ′′(x) ≤ M2

para todo x ∈ [a, b]. La localización del valor de la integral será lo más precisa posible si setoman

M1 = mın{f ′′(x) : x ∈ [a, b]} y M2 = max{f ′′(x) : x ∈ [a, b]}

(si es posible encontrar esos valores extremos, que existen por el teorema del mismo nombre).En cualquier caso se tendría que

−(b− a)3

12 · n2M2 ≤ −(b− a)3

12 · n2f ′′(c) ≤ −(b− a)3

12 · n2M1,

y, por tanto,

h

2

[f(a) + 2

n−1∑i=1

f(xi) + f(b)

]− (b− a)3

12 · n2M2 ≤

∫ b

a

f(x) dx ≤

≤ h

2

[f(a) + 2

n−1∑i=1

f(xi) + f(b)

]− (b− a)3

12 · n2M1.

5.3.2. Error de la fórmula de Simpson

Bajo ciertas hipótesis bastante amplias, el siguiente teorema proporciona una expresióndel error de la fórmula de Simpson.

Teorema 5.2 Sea f una función de clase C4 en [a, b]. Para cierto n ∈ N, sea h = (b−a)/n

y sea a = x0 < x1 < · · · < xn = b la partición de [a, b] de�nida mediante xi = a + ih para

todo i = 0, 1, 2, . . . , n. Consideramos también los puntos medios de los intervalos de dicha

partición, es decir, los puntos xi− 12dados por (5.6). Entonces existe un número c ∈ [a, b] tal

que ∫ b

a

f(x) dx =h

6

[f(a) + 2

n−1∑i=1

f(xi) + 4n∑

i=1

f(xi− 12) + f(b)

]− (b− a)5

2880 · n4f iv(c).

Al ser f iv continua en [a, b], existe M > 0 tal que |f iv(x)| ≤ M para todo x ∈ [a, b].Como en el caso del error de la fórmula de los trapecios, lo ideal será tomar

M = max{|f iv(x)| : x ∈ [a, b]},

siempre que pueda determinarse dicho máximo. Luego el error absoluto cometido por elmétodo de Simpson se puede acotar, ya que∣∣∣∣− (b− a)5

2880 · n4f iv(c)

∣∣∣∣ = (b− a)5

2880 · n4

∣∣f iv(c)∣∣ ≤ M(b− a)5

2880 · n4.

En el último término de la expresión anterior, M , b y a son cantidades �jas. Por tanto,para hacer la cota del error tan pequeña como se desee bastará aumentar el número n desubintervalos del método. El hecho de que esta cota dependa del valor de n4 y de que la cotadel error del método de los trapecios dependiera de n2, indica que el método de Simpson

92

converge a la integral exacta con una rapidez mucho mayor que el método de los trapecios:por ejemplo, al multiplicar el número de subintervalos por 10, la cota del error de los trapeciosse divide por 100, mientras que la cota del error de Simpson se divide por 10000.

Así, si se desea que el error absoluto en el método de Simpson sea menor que ciertacantidad positiva δ, bastará imponer que

M(b− a)5

2880 · n4< δ,

es decir, n deberá satisfacer que

n >4

√M(b− a)5

2880 · δ=

b− a

2

4

√M(b− a)

180 · δ. (5.9)

Veamos un ejemplo de aplicación de esto último.

Ejemplo 5.2 Se desea calcular aproximadamente la integral del ejemplo 5.1 con el mismogrado de aproximación que entonces, pero utilizando la fórmula de Simpson. Para lo cualse determina el número n de subintervalos que asegure que el error no superará dicha cota.Luego n debe satisfacer (5.9) para δ = 10−3. Para determinar la constante M en (5.9) secalcula la cuarta derivada de la función a integrar, que es

f iv(x) = −12 cosx2 + 16x4 cos x2 + 48x2 sen x2.

Puede comprobarse �aunque es un tanto complicado de demostrar� que f v(0) = 0 y quef v(x) > 0 para todo x ∈ (0, 8/9], por lo que f iv es estrictamente creciente en [0, 8/9]. Comof iv(0) = −12 y f iv(8/9) ≃ 25,5286, resulta que

max{|f iv(x)| : x ∈ [0, 8/9]} ≃ 25,5286.

Luego podemos tomar M = 25,5286 y, por tanto,

n >b− a

2

4

√M(b− a)

180 · δ≃ 4

9

4

√25,5286 · (8/9)

180 · 10−3=

4

27

√2 · 25528,6

5≃ 1,4893.

Luego si n ≥ 2 la integral se habrá aproximado con un error inferior a 10−3. 2

En ocasiones se podrá obtener una aproximación algo mejor de la integral si se conocenuna cota inferior M1 y otra superior M2 de f iv(x) en el intervalo [a, b]: M1 ≤ f iv(x) ≤ M2

para todo x ∈ [a, b]. La localización del valor de la integral será lo más precisa posible si setoman

M1 = mın{f iv(x) : x ∈ [a, b]} y M2 = max{f iv(x) : x ∈ [a, b]}(si es posible encontrar estos valores extremos). En cualquier caso se tendría que

− (b− a)5

2880 · n4M2 ≤ − (b− a)5

2880 · n4f iv(c) ≤ − (b− a)5

2880 · n4M1,

y, por tanto,

h

6

[f(a) + 2

n−1∑i=1

f(xi) + 4n∑

i=1

f(xi− 12) + f(b)

]− (b− a)5

2880 · n4M2 ≤

∫ b

a

f(x) dx ≤

≤ h

6

[f(a) + 2

n−1∑i=1

f(xi) + 4n∑

i=1

f(xi− 12) + f(b)

]− (b− a)5

2880 · n4M1.

93

5.4. Ejercicios

1. a) Aproxime mediante el método de los trapecios la integral∫ 2

−1

ex

x2 + 1dx

de modo que solo necesite evaluar el integrando en cuatro puntos.

b) Obtenga una cota del error de aproximación.

2. Se considera la integral∫ 3

1

1

xdx.

a) Aproxime dicha integral mediante el método de los trapecios usando diez subin-tervalos. Compare el resultado con la integral exacta.

b) Obtenga una cota superior y una cota inferior del error cometido por esa aproxi-mación. Compare esas cotas con valor exacto del error.

c) Acote también dicha integral, tanto inferior como superiormente. Compruebe quela integral exacta está entre esas cotas.

d) Repita los apartados anteriores para el método de Simpson, pero usando solocinco subintervalos.

e) ¾Cuántas evaluaciones se necesitarían para asegurar que el error absoluto cometidopor el método de los trapecios fuera inferior a 5 · 10−5? Responda a la mismapregunta para el método de Simpson.

3. Se considera la integral I =

∫ 1

0

1

x2 + 1dx.

a) Aproxime I mediante el método de los trapecios, evaluando la función del inte-grando en 11 puntos.

b) Utilizando la fórmula del error del método de los trapecios, obtenga una cotasuperior y una cota inferior de I.

c) ¾En cuántos puntos sería necesario evaluar el integrando para asegurar que elerror absoluto de la aproximación de I dada por el método de los trapecios fueramenor que 10−4? Obtenga también el número de evaluaciones necesarias para quedicho error sea menor que 10−5 y 10−6.

d) Compare el valor exacto de I con la aproximación calculada en el primer apartadode este ejercicio.

e) Obtenga el error absoluto y el error relativo de esa aproximación.

f ) Como aplicación, encuentre un algoritmo que proporcione aproximaciones delnúmero π tan exactas como se desee.

g) Repita todos los apartados anteriores usando el método de Simpson.

4. Aproxime mediante la fórmula de Simpson la integral de f(x) = e−x2en el intervalo

[2, 4] con un error menor que 2,5 · 10−5.

5. Aproxime la integral de f(x) = cos2 x en el intervalo [0, 1] con un error menor que 10−3

utilizando: (a) el método de Simpson; (b) el método de los trapecios.

94

6. Determine el número de subintervalos que se necesitan a priori, en el método de lostrapecios y en el de Simpson, para aproximar con un error menor que 10−k (k ∈ N) elvalor de la integral ∫ 4

3

esen2 xdx.

7. Use la fórmula de los trapecios para calcular una aproximación de la integral∫ 0,5

0

ecosx dx

con un error absoluto menor que 0,0025. Utilice los mismos puntos usados en el métodoanterior para aproximar la integral mediante la fórmula de Simpson, y obtenga unacota del error cometido.

8. Sea f una función de la que tan solo se conocen los siguientes valores: f(0) = 0,f(1/5) = 1, f(2/5) = 5/2, f(3/5) = 9/2, f(4/5) = 7, f(1) = 10 y f(6/5) = 14.Utilice la fórmula de los trapecios y la fórmula de Simpson para aproximar el valor dela integral de la función f en el intervalo de extremos 0 y 6/5.

9. De una función f se conocen los siguientes puntos de su grá�ca: (1,8, 3,12014),(2, 4,42569), (2,2, 6,04241), (2,4, 8,03014) y (2,6, 10,46675).

a) Utilizando esos datos, obtenga aproximaciones de la integral de f en el intervalo[1,8, 2,6] mediante las fórmulas de los trapecios y de Simpson. ¾Hay algún datoque permita decidir cuál de las dos aproximaciones es la mejor?

b) [Mathematica] Construya el polinomio y el spline cúbico natural que interpolandichos puntos. Úselos para obtener otras dos aproximaciones de dicha integral.

10. Para la realización de un proyecto de ingeniería hidráulica se necesita conocer el área deuna sección de una laguna. En la �gura 5.2 se representan las mediciones de profundidadobtenidas para esa sección. Aproxime su área utilizando la regla de Simpson compuesta.

1.8 2

4 4

6

4 3.63.4 2.8

profundidad

distancia a la orilla izquierda (m.)0 10 20

Figura 5.2: Profundidad de la sección de una laguna

11. La pared de una presa es plana, su base y su borde superior son segmentos horizontalesparalelos y su altura es de 20 m. Se ha medido la anchura de esa pared a nivel del sueloy desde ahí a intervalos de 5 m. de altitud, obteniéndose desde abajo hacia arriba lossiguientes valores: 9, 15, 20, 27 y 30 m. Calcule de forma aproximada el área de lapared de dicha presa mediante la fórmula de Simpson. Indicación: puede suponer para

95

el planteamiento del problema �el resultado en cualquier otro caso sería el mismo�que uno de los bordes laterales de la pared es recto y perpendicular a su base.

12. Dada la función f(x) = sen(πx/4) + 2 sen(πx/2) + 3 cos(πx) − 2 cos(πx/4), aproximesu integral en el intervalo [0, 8] mediante los métodos de los trapecios y de Simpsonde modo que solamente se necesite evaluarla en cinco puntos. Obtenga el polinomiop y el spline cúbico natural s que interpola a f en esos puntos, y aproxime de nuevola integral de f utilizando p y s. Compare las cuatro aproximaciones obtenidas con laintegral exacta.

13. Sea f una función de clase C4 en el intervalo [3, 10]. Sabemos que entonces existe ciertaconstante M > 0 tal que |f iv(x)| ≤ M para todo x ∈ [3, 10].

a) Si se aproxima el valor de la integral de f en dicho intervalo mediante la fórmulade Simpson compuesta, dividiendo el intervalo [3, 10] en n subintervalos iguales,obtenga una cota del error absoluto de la aproximación realizada (en función deM y n).

b) Suponga que se ha tomado cierto valor de n (n = n0), para el que resulta que lacota del error obtenida en el apartado anterior es demasiado grande. Si se necesitaque la cota sea 10000 veces menor, ¾qué se debería hacer?

c) Volviendo al principio del problema, suponga que el error absoluto de aproxi-mación de la integral �al aplicar la fórmula de Simpson compuesta� debe sermenor que 10−8. Obtenga entonces una expresión (en función de M) del númeromínimo de subintervalos de [3, 10] a los que habría que aplicar el método.

14. Se considera una variable aleatoria X que sigue una distribución normal con mediaµ = 0 y desviación típica σ. Como es sabido, su función de densidad es

f(x) =e−x2/2σ2

σ√2π

,

y la probabilidad de que la variable X tome un valor en el intervalo [a, b] es

P [a ≤ X ≤ b] =

∫ b

a

f(x) dx.

Aproxime el valor de las siguientes probabilidades con un error máximo de ±10−5:

(a) P [|X| ≤ σ]; (b) P [|X| ≤ 2σ]; (c) P [|X| ≤ 3σ].

15. Se considera la función g(x) =20

2− x+ x2 + 2x3.

a) Obtenga el polinomio que interpola a g(x) en los puntos x0 = −1, x1 = 0, x2 = 1

y x3 = 2.

b) Utilice ese polinomio de interpolación para dar una aproximación de la integralde g(x) en el intervalo [−1, 2].

c) Utilice también la fórmula compuesta de los trapecios que utiliza los cuatro puntosdados para aproximar esa integral.

96

d) Teniendo en cuenta que −35 ≤ g′′(x) ≤ 95 para todo x ∈ [−1, 2], proporcione unintervalo en el que se encuentre el valor exacto de dicha integral.

e) Suponiendo que dicho valor exacto es 20,0551709, obtenga el error relativo de lamejor de las dos aproximaciones obtenidas en los apartados anteriores.

16. [Mathematica] Se considera la integral∫ 2

−1

ex

x2 + 1dx.

a) Utilizando la instrucción directa de integración numérica de Mathematica, obten-ga una aproximación de esa integral, mostrando al menos diez decimales de laaproximación.

b) Aplique el método de los trapecios para aproximar la integral anterior utilizandon intervalos. Si la aproximación obtenida se denota It(n), obtenga en particularlas aproximaciones It(3), It(30) e It(100).

c) Proporcione una expresión del error cometido por It(3) al aproximar la integral(en función de cierto valor desconocido c ∈ [−1, 2]). Si Et(3, c) denota ese error,represente la función de g(c) = Et(3, c) en el intervalo [−1, 2], y deduzca de ahíuna cota inferior y una cota superior del error de la aproximación It(3).

d) Obtenga el valor de la integral con un error absoluto máximo de 10−6.

17. [Mathematica] Se considera la integral∫ 6

5

√x3 − x− 10

ln(x− 1)dx.

a) Aplique la fórmula de Simpson para aproximar esta integral de manera que lafunción del integrando tenga que evaluarse exactamente en 71 puntos.

b) Acote el error que se comete al aproximar la integral por dicha fórmula.

c) Independientemente de lo anterior, aproxime la integral con un error menor que10−12.

18. [Mathematica] Aplique las fórmulas de los trapecios y de Simpson para aproximar laintegral ∫ 2

1

ln(1− 0,99 cos2

πx

2

)dx,

de manera que la función del integrando tenga que evaluarse exactamente en 101 puntosen ambos casos. Para cada método, acote el error de la aproximación y determine unintervalo en el que se encuentre el valor exacto de la integral. Compare los resultadosanteriores con una aproximación bastante exacta que obtenga mediante la instrucción�NIntegrate�.

19. [Mathematica] Aplique las fórmulas de los trapecios y de Simpson para aproximar laintegral ∫ 2

0

8

x2 + 4dx,

de manera que la función del integrando tenga que evaluarse exactamente en 201 puntosen ambos casos. Para cada método, acote el error de la aproximación y determine unintervalo en el que se encuentre el valor exacto de la integral. Compare los resultadosanteriores con dicho valor exacto.

97

20. a) Dada la función

f(x) =1

x3 + x+ 1,

obtenga la fórmula compuesta de los trapecios �incluyendo la expresión delerror� para la integral ∫ 1

0

f(x) dx,

de modo que solamente se evalúe la función f en seis puntos.b) Sabiendo que 0,25 < f ′′(x) ≤ 2 para todo x ∈ [0, 1], obtenga un intervalo en el

que se encuentre el valor de dicha integral.

21. Se desea aproximar la integral∫ √

3

−√

33

f(x) dx, donde f(x) = arc tg x.

a) Evalúe el integrando en cinco puntos para obtener sendas aproximaciones de laintegral mediante los métodos de los trapecios y de Simpson. Indicación: utiliceque f(−

√33) = −π

6, f(0) = 0, f(

√33) = π

6, f(2

√3

3) = 0,857072 y f(

√3) = π

3.

b) Obtenga una cota del error absoluto de la aproximación obtenida mediante elmétodo de los trapecios.

c) Se ha obtenido que |f iv(x)| < 5 para todo x ∈[−

√33,√3]. Obtenga el número de

puntos en los que se debe evaluar la función f para que se asegure que el métodode Simpson produce una aproximación de la integral cuyo error absoluto es menorque 10−6.

22. Evalúe el integrando de la siguiente integral en siete puntos para obtener sendas aproxi-maciones de la integral mediante los métodos de los trapecios y de Simpson:∫ π

0

√1− sen2 x

2dx.

23. Sea f una función tal que f(0) = 1, f(π/4) = −1, f(π/2) = 2, f(3π/4) = −2,f(π) = 3, f(5π/4) = −3 y f(3π/2) = 4. Utilice la fórmula de los trapecios y la fórmulade Simpson para aproximar el valor de la integral de f en el intervalo [0, 3π/2].

24. La velocidad de un automóvil a los 3, 6, 9 y 12 segundos de ponerse en marcha esrespectivamente de 6, 18, 27 y 33 m/seg.

a) Obtenga el polinomio que interpola los cinco datos dados y estime la velocidaddel automóvil a los 10 segundos desde que se puso en movimiento.

b) Obtenga una aproximación del espacio recorrido por el automóvil en esos 12 se-gundos utilizando: (1) una fórmula de integración numérica de tipo interpolatoriodeducida del apartado a); (2) la fórmula de Simpson.

Algunas cuestiones teóricas

1. ¾Qué signi�ca que una fórmula de integración numérica sea de tipo interpolatorio?

2. En las fórmulas compuestas de los trapecios y de Simpson, si se multiplica por dos elnúmero n de subintervalos utilizados en la fórmula, ¾en qué medida disminuye la cotade error en cada una de esas fórmulas compuestas?

98

Bibliografía

[1] Atkinson, K. E. Elementary Numerical Analysis. Wiley, New York, 1993.

[2] Burden, R. L.; Faires, J. D. Análisis Numérico, 7aedición. Thomson, México, 2002

[3] Ciaurri, O.; Varona, J. L. ¾Podemos �arnos de los cálculos efectuados con ordenador?La Gaceta de la Real Sociedad Matemática Española, 9 (2006), no. 2, 483-514. Tambiénse encuentra en: Boletín de la Sociedad Española de Matemática Aplicada, 37 (2006),93-121.

[4] Chapra, Steven C.; Canale, Raymond P.Métodos numéricos para ingenieros, 4aedición.McGraw-Hill, México, 2005.

[5] Cordero, Alicia; Hueso, José L.; Martínez, Eulalia; Torregrosa, Juan R. Problemas

resueltos de Métodos Numéricos. Thomson, Madrid, 2004.

[6] Faires, J. D.; Burden, R. Métodos numéricos, 3aedición. Thomson, Madrid, 2004.

[7] Gasca, M. Cálculo numérico, 6aedición. Universidad Nacional de Educación a Distan-cia, Madrid, 1988.

[8] Gerald C. F.; Wheatley P. O. Análisis Numérico con Aplicaciones, 6aedición. Pearson,México, 2000.

[9] Kincaid, D.; Cheney, W. Análisis Numérico. Addison-Wesley Iberoamericana, Wilm-ington, 1994.

[10] Ramírez, V.; González, P.; Pasadas, M.; Barrera, D. Matemáticas con Mathematica:Introducción y primeras aplicaciones. Proyecto Sur de Ediciones, Granada, 1996.

[11] Ramírez, V.; Barrera, D.; Pasadas, M.; González, P. Cálculo Numérico con Mathemat-ica. Ariel, Barcelona, 2001.

[12] Rodríguez Lallena, J. A.; López Artés, P.; Martínez González, P. Fundamentos

matemáticos de la Ingeniería Agronómica. Autoedición, Almería, 2002.

[13] Scheid, F.; Di Costanzo, R. E. Métodos numéricos, 2aedición. McGraw Hill, México,1991.

[14] Trott, M. The Mathematica GuideBook for Graphics, Springer, New York, 2004.

[15] Trott, M. The Mathematica GuideBook for Numerics, Springer, New York, 2005.

99

[16] Trott, M. The Mathematica GuideBook for Programming, Springer, New York, 2004.

[17] Trott, M. The Mathematica GuideBook for Symbolics, Springer, New York, 2005.

[18] S. Wolfram, The Mathematica book, 3aedición. Wolfram Media/Cambridge UniversityPress, New York, 1996.

[19] S. Wolfram, Mathematica 3.0 Standard Add-On Packages, Wolfram Media/CambridgeUniversity Press, 1996.

100