Notas de Aula Cálculo Numérico

Notas de Aula

Cálculo Numérico Percy Antonio Ticona Centeno Universidad Nacional de San Agustín Departamento Académico de Matemáticas y Estadística Junio de 2009 Arequipa - Perú

Introducción La matemática está comprendida de varias partes, cada una de ellas tiene importancia propia, pero algunas también son fundamentales en diferentes disciplinas. Muchos problemas de la vida real pueden ser representados por formulaciones matemáticas, las cuales son denominadas modelos matemáticos. Usando argumentos matemáticos teóricos, algunas veces es posible garantizar la existencia de soluciones para esos modelos matemáticos, pero encontrar manualmente esas soluciones puede resultar extremadamente difícil y a veces imposible. Estudiar métodos numéricos desde un punto de vista general nos permitirá analizar mecanismos de cálculo capaces de otorgar soluciones, o aproximaciones a las soluciones, allí donde las herramientas teóricas fracasaban. Estos mecanismos numéricos de cálculo deben caminar de la mano con el computador, pues en su mayoría requieren de muchos pasos y frecuentemente están orientados a la resolución de problemas de grandes dimensiones. La importancia del estudio de métodos numéricos es indiscutible, pues gran parte de las investigaciones en ciencias, ingenierías, economía, etc., recurren a técnicas numéricas para la obtención de resultados.

iii

Capítulo 1 Nociones Básicas Sobre Errores En este primer capítulo aún no estudiaremos los métodos numéricos, sino que veremos algunos conceptos básicos y señalaremos algunos factores que intervendrán en la resolución de problemas mediante el computador. Existen diversas fases cuando intentamos resolver un problema mediante métodos numéricos, la siguiente figura esquematiza ese procedimiento.

Puede suceder que los resultados finales obtenidos no sean justamente los esperados, aunque todas las fases hayan sido ejecutadas correctamente, los motivos pueden ser varios y los estudiaremos a continuación.

1.1.

Factores que intervienen en la resolución numérica de problemas

La mayoría de los métodos numéricos que veremos aquí tienen un carácter iterativo, esto significa que iniciarán con una estimativa inicial solución, © (0)de la ª (0) (1) digamos x , para luego contruir una sucesión de valores x , x , x(2) , ... 1

CAPÍTULO 1. NOCIONES BÁSICAS SOBRE ERRORES

2

de modo que l´ımk→∞ x(k) = x∗ , donde x∗ es la solución. Por lo general, x(k+1) será calculado a partir de x(k) mediante un procedimiento debidamente fundamentado. Frecuentemente, este procedimiento requerirá un gran número de cálculos, por lo que será necesario el auxilio de un computador. Lamentablemente, los computadores presentan algunos inconvenientes que, si no se controlan, pueden ocasionar respuestas catastróficas. Ejemplo 1.1 Considere el trecho de un programa en MatLab, tal como se muestra en la figura 1.1. Observe que en teoría, debería cobrarse $1000000.

Figura 1.1: Un error numérico grave cometido por el computador. Sin embargo, debido al error cometido por el computador, se terminaría pagando $2000000. (Vea también el ejemplo 1.3) Laboratorio 1.1 Haga un programa en algún lenguaje de programación que usted conozca, de modo que en la práctica corrobore el Ejemplo 1.1. Al resolver un problema por métodos numéricos, los resultados obtenidos pueden depender de: 1. La precisión de los datos de entrada 2. La forma cómo éstos son representados en el computador 3. Las operaciones numéricas efectuadas Cada uno de estos temas serán explicados a continuación.

CAPÍTULO 1. NOCIONES BÁSICAS SOBRE ERRORES

1.1.1.

3

Precisión de los datos de entrada

Cuando ingresamos los datos de un problema, ya sea para calcular en un papel o trabajar en el computador, ellos contienen una imprecisión inherente, quiere decir que no hay cómo evitarlos. El siguiente ejemplo aclara esta afirmación. Ejemplo 1.2 Sabemos que para calcular el área de un círculo, tenemos que ingresar numéricamente el radio r y el valor de π. El valor de r quizá pueda ser conocido exactamente (r = 2), pero apenas podemos conocer una aproximación de π con un número finito de dígitos. Así, aproximando el valor de π por 3,14, el área del círculo será: 3,14 × (2)2 = 12,56 m2 Si consideramos 3,1416, entonces el área del círculo estará dada por: 3,1416 × (2)2 = 12,5664 m2 Pero si consideramos 3,141592654, entonces el área será 3,141592654 × (2)2 = 12,566370616 m2 Claramente, las imprecisiones de los datos de entrada ocasionan imprecisiones en los resultados. En el ejemplo anterior vimos que el mejor resultado se obtuvo en el último caso. Pero, cuando usamos un computador, ¿cuántos dígitos decimales reconoce éste? El siguiente ejemplo intentará aclarar esta situación. Ejemplo 1.3 Usando MatLab 6.0 en una PC de 32-bit con sistema operativo Windows XP, hicimos la siguiente operación: 0,00000000000001 + 1 y el resultado obtenido fue 1,0000000000001, lo cual es satisfactorio. Luego, hicimos: 0,000000000000001 + 1 y obtuvimos como respuesta 1. ¿Qué sucedió? El ejemplo anterior nos hace ver que en el primer caso, cuando se usan los 14 dígitos decimales a la derecha del punto, el sistema de cómputo no comete error. Mientras que, si se usan 15 dígitos, el sistema nos otorga una respuesta errada. La razón se debe a que todo computador trabaja con un número finito y bien reducido de dígitos, en nuestro caso 14 a la derecha del punto, si el número de dígitos sobrepasa lo esperado, el sistema lo trunca o redondea, dependiendo del sistema utilizado.

CAPÍTULO 1. NOCIONES BÁSICAS SOBRE ERRORES

4

Laboratorio 1.2 Utilice un paquete o lenguaje de programación en un computador, para corroborar los resultados obtenidos en el Ejemplo 1.3.

1.1.2.

Representación de los datos en el computador

Un número con representación decimal finita puede tener una representación infinita en el sistema binario. Como un computador trabaja con el sistema binario y con una cantidad fija de dígitos, necesariamente trabajará con una aproximación, por lo tanto no se obtendrán resultados finales exactos. El sistema con el que trabajamos comunmente es el decimal, un número x en este sistema lo representaremos algunas veces, cuando se preste a confusión, por (x)10 . Por otro lado, el sistema con el que trabaja un computador hoy en día es el binario, un número y en este sistema será representado por (y)2 . Así por ejemplo, (5)10 y (101)2 representan el número 5 en el sistema decimal y binario, respectivamente. El siguiente ejemplo muestra cómo esta representación aparentemente inofensiva, puede generar terribles errores cuando se trabaja en un computador. Ejemplo 1.4 Considere la siguiente sumatoria: S=

30000 X

ai

i=1

Para ai = 0,5, el resultado exacto debería ser S = 0,5

30000 X i=1

1 = 0,5 × 30000 = 15000

Después de implementar un pequeño programa en computador, el resultado fue también 15000. Claramente, no hay por qué preocuparse en este caso, los resultados son los mismos. Pero, para ai = (0,11)10 , el resultado exacto debería ser S = (0,11)10 × 30000 = 3300 Sin embargo, el resultado obtenido por el computador fue S = 3300,00000000063 ¿Cómo explicar la diferencia de resultados en este caso?

CAPÍTULO 1. NOCIONES BÁSICAS SOBRE ERRORES

5

Esto se debe a que el número (0,11)10 , cuya representación en el sistema decimal es finita, tiene una representación binaria infinita: ¢ ¡ 0,000111000010100011110101110000101000111101 2

Si el computador trabajara con 14 dígitos después del punto, el número debería ser cortado o redondeado, lo cual representa ya un error. Todos los cálculos subsiguientes serán afectados por este hecho.

Laboratorio 1.3 Utilizando algún lenguaje de programación, haga un programa para ejecutar lo tratado en el ejemplo 1.4. Por lo general, el error ocurrido depende de la representación de los números en la máquina utilizada. La representación de un número depende de la base elegida o disponible en la máquina en uso, y, del número máximo de dígitos usados en su representación. Cualquier cálculo que envuelva números que no pueden ser representados a través de un número finito de dígitos, no otorgará como resultado un valor exacto. Cuanto mayor sea el número de dígitos utilizados, mayor será la precisión obtenida. Como vimos en el ejemplo 1.4, un número puede tener representación finita con respecto a una base, pero una representación infinita en otra base. La base decimal es la que empleamos generalmente, pero antiguamente fueron empleadas otras bases, como la base 12 y la base 60. Un computador opera normalmente en el sistema binario. Observe lo que pasa cuando un usuario interactúa con el computador: Los datos de entrada son enviados al computador por el usuario en el sistema decimal, esa información es convertida al sistema binario por el computador, y, todas las operaciones son realizadas en ese sistema. Los resultados finales serán convertidos para el sistema decimal y, finalmente, serán transmitidos hacia el usuario. Todo este proceso es una fuente de errores que afecta el resultado final de los cálculos.

1.1.3.

Las operaciones numéricas efectuadas

Pero errores no sólo ocurren el la imprecisión de los datos de entrada y su representación binaria. Errores ocurren también en las operaciones numéricas efectuadas por un sistema de cómputo (binario).

CAPÍTULO 1. NOCIONES BÁSICAS SOBRE ERRORES

6

Conversión de Números en los Sistemas Decimal y Binario Cualquier número entero en la base β, de la forma (aj aj−1 , ...a2 a1 a0 )β , donde 0 ≤ ak ≤ β − 1 y k = 0, ..., j, puede ser escrito en la forma polinomial aj β j + aj−1 β j−1 + ... + a2 β 2 + a1 β 1 + a0 β 0 Mediante esa representación, podemos convertir fácilmente un número entero representado en el sistema binario para el sistema decimal, e inversamente. Por ejemplo: (10111)2 puede ser representado por 1 × 24 + 0 × 23 + 1 × 22 + 1 × 21 + 1 × 20 Reordenando y resaltando la base 10 (10111)2 = 1 × 24 + 0 × 23 + 1 × 22 + 1 × 21 + 1 × 20 ¡ ¢ = 2 23 + 2 + 3 = 2 × 101 + 3 × 100 = (23)10 Laboratorio 1.4 En algún lenguaje de programación, haga un programa tal que, dado un número entero binario, retorne su equivalente decimal. E inversamente, dado un entero decimal, otorgue su equivalente binario. ¿Cómo Convertir un Número Fraccionario de Representación Decimal a Binario? Consideremos ahora la conversión de un número fraccionario de base 10 para la base 2. Por ejemplo, r = 1,25, s = 0,666..., t = 0,414213562..., etc. Notemos que r tiene una representación finita, pero s y t tienen representaciones infinitas. En términos generales, dado un número entre 0 y 1 en el sistema decimal, ¿cómo obtener su representación binaria? Considerando el número decimal fraccionario r = 0,125, existen dígitos binarios d1 , d2 , ..., dj , ... tal que (0.d1 d2 ...dj ...)2 será su representación binaria en la base 2. Así, (0,125)10 = d1 × 2−1 + d2 × 2−2 + ... + dj × 2−j + ... Multiplicando cada término de la expresión (1.1) por 2, tendremos 2 × 0,125 = d1 + d2 × 2−1 + d3 × 2−2 + ... + dj × 2−j+1 + ...

(1.1)

CAPÍTULO 1. NOCIONES BÁSICAS SOBRE ERRORES

7

Por tanto, d1 representa la parte entera de 2×0,125, que es igual a 0, mientras que d2 × 2−1 + d3 × 2−2 + ... + dj × 2−j+1 + ... representa la parte fraccionaria de 2 × 0,125, que es 0,250. Aplicando ahora el mismo procedimiento para 0,250, tendremos 0,250 = d2 × 2−1 + d3 × 2−2 + ... + dj × 2−j+1 + ... 2 × 0,250 = 0,5 = d2 + d3 × 2−1 + d4 × 2−2 + ... + dj × 2−j+2 + ... de donde d2 = 0. Repitiendo el procedimiento para 0,5 tenemos 0,5 = d3 × 2−1 + d4 × 2−2 + ... + dj × 2−j+2 + ... 2 × 0,5 = 1 = d3 + d4 × 2−1 + ... + dj × 2−j+3 + ... de donde d3 = 1. Como la parte fraccionaria de 2 × 0,5 es cero, el proceso de conversión termina. En resumen tenemos: d1 = 0, d2 = 0 y d3 = 1. Por lo tanto, el número (0,125)10 tiene representación finita en la base 2: (0,125)10 = (0,001)2 Laboratorio 1.5 Usando los procedimientos anteriores para convertir números fraccionarios decimales, a binarios, haga un programa usando algún lenguaje de programación y verifique que: 1. El número (0,5)10 tiene una representación binaria finita (0,1)2 2. El número (0,11)10 tiene una representación binaria infinita ¢ ¡ 0,000111000010100011110101110000101000111101 2 3. Verifique cuántos dígitos el computador está considerando.

El hecho de que un número no tenga representación finita en el sistema binario, puede ocasionar la ocurrencia de errores aparentemente inexplicables en los cálculos efectuados en sistemas computacionales binarios. Un computador que opera en el sistema binario, necesariamente tendrá que almacenar una aproximación para (0,11)10 , debido a que sólo posee una cantidad fija de posiciones para guardar los dígitos de la mantisa de un número. Al ser esta aproximación usada para realizar los cálculos, no puede esperarse un resultado exacto.

CAPÍTULO 1. NOCIONES BÁSICAS SOBRE ERRORES

8

Laboratorio 1.6 (Precisión de una máquina) La precisión de la máquina se define como el menor número ε > 0 en aritmética de punto flotante, tal que (1 + ε) > 1. Este número depende totalmente del sistema de representación de la máquina: base numérica, total de dígitos en la mantisa, de la forma cómo son realizadas las operaciones y del compilador utilizado. Es importante conocer la precisión de la máquina porque en varios algoritmos se requiere ingresar como dato de entrada un valor positivo, próximo de cero, para que sea usado como criterio de comparación para la detención del algoritmo. El siguiente algoritmo determina dicha precisión: Paso 1: A ← 1, S ← 2 Paso 2: Mientras S > 1, hacer A 2 S ← 1+A

A ←

Paso 3: Hacer prec = 2A e imprimir prec. 1. Haga un programa en algún lenguaje de programación que ejecute el algoritmo anterior. 2. Discuta su significado práctico. Aritmética de Punto Flotante Cualquier computador o calculadora representa un número en un sistema denominado aritmética de punto flotante. En este sistema, el número r será representado en la forma: ± (.d1 d2 ...dt ) × β e donde β t e

: La base en que la máquina opera : El número de dígitos en la mantisa, 0 ≤ dj ≤ β − 1, j = 1, ..., t, d1 6= 0 : El exponente en el intervalo [−u; u], el valor de u depende de la máquina con que se esté trabajando

En una computadora, sólo una pequeña cantidad de números son representados exactamente, por lo general, la representación será realizada por medio de truncamiento o redondeo.

CAPÍTULO 1. NOCIONES BÁSICAS SOBRE ERRORES

9

Ejemplo 1.5 Considere una máquina que opera en el sistema β = 10, t = 3, e ∈ [−5; +5]. Los números no nulos representados en este sistema serán de la forma ± (.d1 d2 d3 ) × 10e ,

0 ≤ dj ≤ 9,

d1 6= 0,

e ∈ [−5; +5]

El menor número, no nulo y en valor absoluto, expresado en esta máquina será: m = 0,100 × 10−5 Mientras que el mayor número es:

M = 0,999 × 10+5 Ahora, en esta misma máquina. consideremos el subconjunto de números reales caracterizados por: G = {x ∈ R : m ≤ |x| ≤ M} Pueden ocurrir varias cosas: x ∈ G Por ejemplo, si x = 235,89 = 0,23589 × 10+3 . Este número posee 5 dígitos en la mantisa. Debido a que t = 3, este número no será considerado de forma exacta en esta máquina. Si la máquina usa truncamiento, entonces el número será representado como 0,235 × 10+3 . Pero si la máquina usa redondeo, entonces el número será representado por 0,236 × 10+3 . |x| < m Por ejemplo, si x = ±0,345×10−7 . Este número no puede ser representado en esta máquina porque el exponente e es menor que −5. La máquina en estas condiciones retorna una advertencia de underflow. |x| > M Por ejemplo, x = ±0,875×109 . En este caso, el exponente es mayor que 5 y la máquina no lo puede representar, advierte la ocurrencia de overflow. Algunos lenguajes de computador permiten que las variables sean declaradas en doble precisión. En este caso, tal variable será representada en el sistema de aritmética de punto flotante de la máquina, pero con aproximadamente el doble de dígitos disponibles en la mantisa. Debemos resaltar que en estas condiciones, el tiempo de ejecución y requerimientos de memoria aumentan considerablemente. La adición en aritmética de punto flotante requiere el alineamiento de los puntos decimales de los dos números. Para eso, la mantisa de menor exponente debe ser desplazada para la derecha.

CAPÍTULO 1. NOCIONES BÁSICAS SOBRE ERRORES

10

Ejemplo 1.6 Sumar x = 0,937 × 104 y y = 0,1272 × 102 Alineando los puntos decimales, tenemos x = 0,937 × 104 y y = 0,001272 × 104 Entonces, x + y = (0,937 + 0,001272) × 104 = 0,938272 × 104 El resultado almacenado después del truncamiento será 0,9382×104 . Mientras que después del redondeo será 0,9383 × 104 . El cero puede representarse con una mantisa nula y cualquier exponente. Por lo general, se utiliza el menor exponente posible de la máquina. Caso contrario, si se usa cualquier exponente para denotar el cero, se pueden perder dígitos significativos, tal como muestra el siguiente ejemplo. Ejemplo 1.7 Supongamos que tenemos una máquina que opera con base 10 y 4 dígitos en la mantisa. Si denotáramos al cero por 0,0000×104 , al sumarlo al número y = 0,3134 × 102 : 0,0000 × 104 + 0,3134 × 102 = 0,0000 × 104 + 0,003134 × 104 = (0,0000 + 0,003134) × 104 = 0,003134 × 104 El resultado después del truncamiento sería 0,0031×104 = 0,3100×102 . Esto significa que fueron perdidos 2 dígitos del valor exacto. Ejemplo 1.8 Represente los siguientes números en un sistema de aritmética de punto flotante (con redondeo y con truncamiento) de 3 dígitos, cuando β = 10, m = 10−4 y M = 10+4 : 3,14

10,053

− 238,15

2,71828

0,000007

718235,82

Para el primer caso, con truncamiento nos resulta 0,314 × 10, mientras que con redondeo 0,314 × 10. Para el segundo caso, con trucamiento obtenemos 0,100 × 102 y con redondeo 0,100 × 102 . Y así sucesivamente: −0,238 × 103 0,271 × 10 underflow overflow

−0,238 × 103 0,272 × 10 underflow overflow

CAPÍTULO 1. NOCIONES BÁSICAS SOBRE ERRORES

11

Todo esto nos da una idea de los posibles errores que pueden suceder, ya sea por desconocimiento nuestro o por la limitación del computador, en el proceso de la resolución numérica de problemas. Debemos advertir que aún es posible realizar un análisis más completo del manejo de errores, pero eso lo veremos en otra ocasión.

1.2.

Errores Absolutos y Errores Relativos

Definimos el error absoluto como la diferencia entre el valor exacto de un número x y su valor aproximado x¯: EAx = x − x¯

(1.2)

Frecuentemente, sólo nos interesa la magnitud de este error. Así por ejemplo, si x = 12,60 y x¯ = 12,81, el error absoluto es EAx = 12,60 − 12,81 = −0,21. Mientras que la magnitud de este error es |−0,21| = 0,21. Esta idea puede ser extendida para comparar la proximidad de vectores. Por ejemplo, consideremos los vectores en R3 dados por     6,1 6,2 x =  5,8  y x¯ =  5,9  −11,3 −10,9

Entonces, usando la norma euclídea, la magnitud del error estaría dada ahora por q kx − x¯k = (6,1 − 6,2)2 + (5,8 − 5,9)2 + (−11,3 − (−10,9))2 = 0,4243

No obstante, el error absoluto definido en (1.2) quizá no tenga interés práctico en este caso.

En general, apenas el valor de x¯ es conocido, lo que hace imposible obtener el error absoluto exacto. Lo que se puede hacer en ese caso es obtener una cota superior o una estimativa para el módulo del error absoluto, tal como muestra el siguiente ejemplo. Ejemplo 1.9 Conociéndose que π ∈ h3,14; 3,15i, podemos tomar para x¯ un valor dentro de este intervalo y tendremos que: |EAπ | = |π − x¯| < 0,01 En este caso diremos que el error absoluto de x¯ con respecto a π, en módulo, es menor que 0,01. Más aún, diremos que el número x¯ está representado con presición menor que 0,01.

CAPÍTULO 1. NOCIONES BÁSICAS SOBRE ERRORES

12

Ejemplo 1.10 En la práctica, a veces no es recomendable controlar algunos procesos basados en el error absoluto. Por ejemplo, si usted gana un premio de S/. 10000000, y cuando va a recogerlo le dicen que sólo tienen S/. 9999990, entonces a usted puede que no le importe la diferencia, pues apenas hay un error absoluto de S/. 10. Pero qué pasa si usted ganó S/. 20 de premio, y cuando usted va a recogerlo le dicen que apenas tienen S/. 10, observe que el error absoluto sigue siendo S/. 10, probablemente no le agrade nada esta última situación, pues se trata de la mitad del premio. Para evitar situaciones como la anterior, en la práctica es mejor utilizar otro criterio para medir el error, éste es conocido como error relativo. El error relativo es definido como el valor absoluto dividido por el valor aproximado, es decir: x − x¯ EAx ERx = = x¯ x¯ Frecuentemente, también se suele trabajar con el módulo de este valor. Observe que el error relativo respecto al primer premio de S/.10000000 es 10 ≈ 0,000001 9999990 Mientras que el error relativo respecto al segundo premio es ER10000000 =

10 =1 10 Con esto, digamos que en este caso se midió el error con más justicia. ER20 =

Ejercicio 1.1 Convierta los siguientes números decimales para su forma binaria: 26, 1278 y 0,1217. Ejercicio 1.2 Convierta los siguientes números binarios para su forma decimal: (101101)2 , (0,111111101)2 y (0,1101)2 . Laboratorio 1.7 El siguiente algoritmo calcula de una forma aproximada la raíz n-ésima de un número no negativo a. Ingresan: a, n y ε > 0, donde ε es la precisión deseada (ε = 10−9 ) x=a Mientras |xn − a| > ε (controlando la magnitud del error absoluto) x = x − (xn − a) / (nxn−1 ) , x 6= 0 √ n Retorna x (una aproximación para a) En algún lenguaje de programación, haga un programa para ejecutar este algoritmo. Modifique el programa para que retorne también el número de pasos (iteraciones). ¿Cómo utilizaría el error relativo para controlar el algoritmo?

CAPÍTULO 1. NOCIONES BÁSICAS SOBRE ERRORES

13

Laboratorio 1.8 (Cálculo de ex ) En algún lenguaje de programación, haga un programa para calcular ex mediante la serie de Taylor con n términos. El valor de x y el número de términos de la serie, n, deben ser dados en la entrada de su programa. Para valores negativos de x, el programa debe calcular ex de dos formas: En una de ellas el valor de x es usado directamente en la serie de Taylor y, en la otra forma, el valor usado en la serie será y = −x, y en seguida, se calcula el valor de ex por medio de e1x . 1. Experimente su programa con varios valores de x (x próximo de cero y distante de cero) y, para cada valor de x, experimente el cálculo de la serie con varios valores de n. Analice los resultados obtenidos. 2. (Dificultades con el cálculo del factorial) El cálculo de k! necesario en la serie de Taylor puede ser hecho de modo a evitar la ocurrencia de overflow. Para esto es necesario analizar cuidadosamente el k-ésimo k término, xk! . Intente combinar el cálculo del numerador con el del denominador y realizar divisiones intermedias. Estudie una manera de realizar esta operación de modo que k! no se sobrecargue. 3. Con la modificación del segundo ítem, la serie de Taylor puede ser calculada con los términos que se desee. ¿Cuál sería el criterio para detener su programa e interrumpir el cálculo de la serie?

Capítulo 2 Ceros reales de funciones reales En esta sección vamos a resolver la ecuación representada por f (x) = 0

(2.1)

donde f : [a; b] ⊂ R 7→ R. Resolver tal ecuación significa encontrar ξ ∈ [a; b] de modo que f (ξ) = 0. Algunas de las técnicas para resolver esta ecuación, es decir, encontrar una raíz de la ecuación (2.1) o simplemente un cero de f , requieren de un procedimiento que comprende esencialmente dos fases. En la primera fase un intervalo conteniendo la raíz es obtenido (aislamiento de las raíces). En la segunda fase, se obtiene una aproximación de la raíz deseada (refinamiento).

2.1.

Aislamiento de las Raíces

Si ξ es una raíz de f , el procedimiento de aislamiento de una raíz consiste en obtener un intervalo [a; b] que contenga ξ. Una primera alternativa sería mediante una observación gráfica, probablemente con ayuda de un computador o una calculadora. √ Ejemplo 2.1 Para el aislamiento de una raíz de f (x) = x − 5e−x , procedemos a graficar utilizando MatLab. Observamos que en el intervalo [0; 4] se encuentra una raíz. La figura 2.1 muestra esta situación. Otra alternativa para el aislamiento de una raíz consiste en analizar el cambio de signo de los valores de la función, tal como explica el siguiente ejemplo. 14

CAPÍTULO 2. CEROS REALES DE FUNCIONES REALES

15

Figura 2.1: Aislando gráficamente una raíz con ayuda de MatLab. √ Ejemplo 2.2 Considere nuevamente la función f (x) = x − 5e−x . Determine un intervalo donde se encuentre por lo menos una raíz. Solución: Tenemos que Df = R+ 0 (los reales positivos y el cero), construimos una tabla de valores con el signo de f (x) para determinados valores de x: x 0 1 2 3 ··· f (x) − − + + · · · Analizando la tabla, vemos que f admite por lo menos una raíz en el intervalo [1; 2]. Para ver si esa raíz es única en ese intervalo, podemos analizar el signo de la derivada de f : 1 f 0 (x) = √ + 5e−x > 0, 2 x

∀x > 0

Vemos que f es estrictamente creciente en R+ . Por lo tanto, podemos concluir que f admite una única raíz en Df y ésta se encuentra en el intervalo [1; 2]. Ejercicio 2.1 Sea f : R 7→ R una función continua en R. Si sabemos que una raíz ξ de f está en [a; b]. Qué podría usted decir sobre ξ si: 1. f (a) f (b) < 0 2. f (a) f (b) > 0

CAPÍTULO 2. CEROS REALES DE FUNCIONES REALES

16

3. f (a) f (b) ≤ 0 4. f (a) f (b) ≥ 0 Laboratorio 2.1 Usando un computador y un programa adecuado, grafique las siguientes funciones y determine aquellos intervalos que incluyan alguna raíz: 1. f (x) = 1000 sin(x3 + 1)/ log(5 + x2 ), x ∈ [−2; 1] 2. g(x) = ln (x2 + 1) − 200 (x + 10)3 + 9x2 + 5, x ∈ R 3. h(x) = 0,00037x11 − (x − π)2 + x2 + 5x − 100, x ∈ R

2.2.

Refinamiento: Métodos Iterativos

Una vez que tenemos un intervalo que contenga la raíz, el siguiente paso es construir un mecanismo que nos otorgue aproximaciones razonables a la solución exacta. En esta sección veremos algunos métodos numéricos clásicos los cuales nos otorgarán aproximaciones a una raíz de f . En su mayoría, estos métodos son iterativos, es decir, que inician con una estimativa de la solución inicial y utilizan ésta para encontrar la siguiente, y así por delante, hasta obtener una aproximación a la solución. Cuando se use un método iterativo, debemos considerar un criterio para detener el algoritmo respectivo. En los métodos que buscan una raíz, éstos se repetirán hasta que xk sea próxima a la raíz exacta ξ con precisión ε > 0, esto ocurrirá si: 1. |xk − ξ| < ε, o ¡ ¢ 2. |f xk | < ε

Pero, ¿cómo efectuar el primer ítem si no se conoce ξ? Una forma es reducir el intervalo que contiene a la raíz en cada iteración. Al conseguirse un intervalo [a; b] tal que ξ ∈ [a; b]

y

b−a 0, hacer a = c e ir al paso 5. Caso contrario, hacer b = c e ir al paso 5. 5. Si b − a < ε, elegir x¯ ∈ [a; b] y finalizar el algoritmo. Caso contrario, hacer k = k + 1 y volver al paso 3. Estudio de la Convergencia del Método de Bisección Bajo las hipótesis establecidas, es claro que el método de la bisección construirá una sucesión {xk } que converge a una raíz. Para probar esto analíticamente procedemos del siguiente modo. Supongamos que [a0 ; b0 ] sea el intervalo inicial y ξ la única raíz de f en el interior de ese intervalo. El método de la bisección genera tres sucesiones: La sucesión {ak }k∈N : No decreciente y acotada superiormente por b0 . Luego, existe r ∈ R tal que l´ımk→∞ ak = r. La sucesión {bk }k∈N : No creciente y limitada inferiormente por a0 . Luego, existe s ∈ R tal que l´ımk→∞ bk = s. La sucesión {ck }k∈N : Generada por ck = todo k = 1, 2, ....

ak +bk , 2

donde ak < ck < bk , para

Observe que el tamaño de cada intervalo es la mitad del intervalo anterior. Así, para k = 1, 2, ... b0 − a0 2 b0 − a0 b1 − a1 = = 2 22 b2 − a2 b0 − a0 = = 2 23 .. . b0 − a0 = 2k

b1 − a1 = b2 − a2 b3 − a3 bk − ak

(2.2)

CAPÍTULO 2. CEROS REALES DE FUNCIONES REALES

20

Entonces, (b0 − a0 ) =0 k→∞ k→∞ 2k y {bk }k∈N son convergentes: l´ım (bk − ak ) = l´ım

Como {ak }k∈N

l´ım bk − l´ım ak = 0 =⇒ s − r = 0

k→∞

k→∞

Por lo tanto, s = r. Sea λ = s = r el límite de las dos sucesiones. Dado que para todo k = 1, 2, ... el punto ck ∈ hak ; bk i, entonces l´ım ck = λ

k→∞

Resta probar apenas que λ es un cero de la función f , o sea, f (λ) = 0. Observe que en cada iteración k = 1, 2, ... tenemos que f (ak ) f (bk ) < 0. Entonces, dado que hemos asumido f continua en [a; b]: 0 ≥

l´ım f (ak ) f (bk ) = l´ım f (ak ) l´ım f (bk )

k→∞

k→∞

k→∞

= f ( l´ım ak )f ( l´ım bk ) = f (r) f (s) = f (λ) f (λ) k→∞

k→∞

= (f (λ))2 de donde concluimos que f (λ) = 0. Por tanto, l´ımk→∞ ck = λ es un cero de f , como habíamos asumido que en el intervalo había una única raíz, tenemos que λ = ξ. Estimación del Número de Iteraciones Dada una precisión ε > 0 y un intervalo inicial [a; b], es posible saber, a priori, cuántas iteraciones serán efectuadas por el método de la bisección (algoritmo 2.1) hasta que se cumpla b − a < ε. Vimos en (2.2) que para k = 1, 2, ... bk − ak =

b0 − a0 bk−1 − ak−1 = 2 2k

(2.3)

Observe que el algoritmo 2.1 se detendrá cuando bk −ak < ε, según la ecuación (2.3) esto equivale a encontrar un valor de k de modo que b0 − a0

b0 − a0 ε

Lo cual implica que k > log2 (b0 − a0 ) − log2 (ε)

(2.4)

Por lo tanto, en el algoritmo 2.1, si el número de iteraciones k es por lo menos blog2 (b0 − a0 ) − log2 (ε)c + 1 el intervalo [ak ; bk ] conteniendo a la raíz ξ verifica bk − ak < ε, en estas condiciones cualquier x¯ ∈ [ak ; bk ] debería ser la aproximación a la raíz buscada ξ, pues |¯ x − ξ| ≤ bk − ak < ε. Podemos concluir sobre este algoritmo que el número de iteraciones dependerá de la longitud del intervalo [a; b] y de la precisión deseada ε. Como muestra la ecuación (2.4), ese número de iteraciones no debería ser grande, debido a la presencia del logaritmo. Ejemplo 2.3 Usando el algoritmo 2.1, se desea encontrar una aproximación a un cero de la función definida por f (x) = x log10 x − 1 la cual está en el intervalo [2; 3] y con precisión ε = 10−2 . ¿Cuántas iteraciones debería efectuar el algoritmo? Solución: Según (2.4) vemos que ¥ ¡ ¢¦ ¥ ¡ ¢¦ k ≥ log2 (3 − 2) − log2 10−2 + 1 = log2 1 − log2 10−2 + 1 ¡ ¢¦ ¥ = − log2 10−2 + 1 = 7 Luego, el algoritmo debería detenerse con k = 7 iteraciones.

Conclusión 2.1 (sobre el método de bisección) Si f : R 7→ R es continua en el intervalo [a; b] y f (a) f (b) < 0: El método de bisección genera una sucesión que converge a la raíz de f (x) = 0, esto se consigue mediante las reducciones sucesivas del intervalo de búsqueda hasta una precisión deseada. Cada iteración no requiere cálculos complicados.

CAPÍTULO 2. CEROS REALES DE FUNCIONES REALES

22

No se requieren derivadas. Las hipótesis no son rigurosas. Converge razonablemente rápido (comparado a otros métodos se le considera lento). Laboratorio 2.2 En algún lenguaje de programación de su preferencia, implemente el algoritmo de la bisección, resuelva las siguientes ecuaciones: 1. ln (x2 + 1) = 200 (x + 10)3 − 9x2 − 5, en R 2. 1000 sin(x3 + 1)/ log(5 + x2 ) = 0, en [−2; 1] 3. 0,00037x11 − (x − π)2 + x2 = −5x + 100, en R

2.2.2.

Método de la Posición Falsa

Sea f : R 7→ R continua en [a; b] tal que f (a) f (b) < 0. Suponga que el intervalo [a; b] contiene una única raíz de la ecuación f (x) = 0. Podemos esperar conseguir una raíz aproximada usando las informaciones sobre los valores de f disponibles a cada iteración. Por ejemplo, en la figura 2.4 se aprecia que al ser |f (a)| pequeño en comparación a |f (b)|, podemos sospechar que la raíz se encuentra más cercana al punto a que al punto b. Luego, en cada iteración, en vez de tomar ck como el punto medio, como lo hacía el método de bisección, podemos tomarlo de la siguiente manera: a |f (b)| + b |f (a)| ck = |f (b)| + |f (a)| que en realidad es una media aritmética ponderada entre a y b, con pesos |f (b)| y |f (a)|. Después de unos cálculos, tenemos: ck =

af (b) − bf (a) f (b) − f (a)

Lo que resta del método de la posición falsa es análogo al método de bisección, la parte donde no se encuentra la raíz debería ser desechada y el intervalo debería ser reducido hasta una presición deseada. Lamentablemente, las cosas no son como parecen, pues el valor de |f (a)| puede ser pequeño y sin embargo, la raíz puede estar muy lejos de a, tal como se aprecia en la figura 2.5. Esto muestra que la sospecha puede estar totalmente errada, lo que retrasaría la convergencia del método, de ahí el nombre.

CAPÍTULO 2. CEROS REALES DE FUNCIONES REALES

23

Figura 2.4: Éxito en el método de la posición falsa: |f (a)| pequeño y raíz cercana de a.

Figura 2.5: Fracaso en el método de la posición falsa: |f (a)| pequeño y raíz alejada de a.

CAPÍTULO 2. CEROS REALES DE FUNCIONES REALES

24

Ejercicio 2.2 Haga un algoritmo que resuma el método de la posición falsa. Laboratorio 2.3 En algún lenguaje de programación de su preferencia, implemente el algoritmo de la posición falsa, pruebe con varios ejemplares y compare el número de iteraciones con el método de bisección. Ejercicio 2.3 Investigue sobre la convergencia y el número de iteraciones requeridas por el método de la posición falsa.

2.2.3.

Método del Punto Fijo (MPF)

Sea f : R 7→ R una función. Se dice que λ es un punto fijo de f , si f (λ) = λ. Este concepto es el mismo para una función vectorial de variable vectorial. El Método del Punto Fijo consiste en lo siguiente: 1. Transformar la ecuación f (x) = 0 en una ecuación equivalente: x = θ (x) 2. Dado un punto inicial x0 ∈ R, generar una sucesión {xk } de aproximaciones hacia ξ, la raíz buscada, mediante la relación xk+1 = θ (xk ) Observe que la función θ debe ser una función que cumpla: f (ξ) = 0 si, y sólo si, θ (ξ) = ξ. De este modo, resolver el problema de encontrar una raíz de una ecuación se convierte en un otro problema de hallar un punto fijo de una función. Aunque a simple vista no parezca, más adelante veremos que esta idea trae ciertas ventajas. La función θ con esa característica se denomina función iteración asociada a la ecuación f (x) = 0. Naturalmente, pueden existir muchas funciones iteración asociadas a una sola ecuación. Ejemplo 2.4 Dada la ecuación x2 + 2x − 10 = 0. Algunas candidatas a función iteración son las siguientes: 1. θ (x) = 5 − 2. θ (x) =

10 , x+2

x2 2

x 6= −2

CAPÍTULO 2. CEROS REALES DE FUNCIONES REALES 3. θ (x) =

10 x

− 2,

25

x 6= 0

Definición 2.1 (Forma general de una función iteración) Una función iteración asociada a la ecuación f (x) = 0, está dada de una forma general por: θ (x) = x + A (x) f (x) con la condición que en ξ, punto fijo de θ, se tenga A (ξ) 6= 0. Teorema 2.1 Si θ es una función iteración de la ecuación f (x) = 0, entonces f (ξ) = 0 si, y sólo si, θ (ξ) = ξ. Prueba. ( =⇒ ) Sea ξ tal que f (ξ) = 0. Como θ (ξ) = ξ + A (ξ) f (ξ), claramente tenemos θ (ξ) = ξ. ( ⇐= ) Si θ (ξ) = ξ, entonces ξ + A (ξ) f (ξ) = ξ implica que A (ξ) f (ξ) = 0, esto a su vez implica que f (ξ) = 0, pues A (ξ) 6= 0 por definición. Estudio de la Convergencia del Método del Punto Fijo Dependiendo de la elección de la función iteración θ, el método del punto fijo puede o no convergir a la solución de la ecuación f (x) = 0. El siguiente teorema establece las condiciones suficientes para que esta convergencia suceda. Teorema 2.2 Sea ξ una raíz de la ecuación f (x) = 0, aislada en un intervalo abierto I centrado en ξ. Sea θ una función iteración asociada a esta ecuación. Si 1. θ y θ0 son funciones continuas en I. 2. |θ0 (x) | ≤ M < 1, para todo x ∈ I. 3. x0 ∈ I Entonces, la sucesión {xk } generada por la regla xk+1 = θ (xk ), k = 0, 1, 2, ... converge hacia ξ. Prueba. La prueba consta de dos partes: 1. Si x0 ∈ I, entonces xk ∈ I, para todo k = 0, 1, 2, ... 2. l´ımk→∞ xk = ξ

CAPÍTULO 2. CEROS REALES DE FUNCIONES REALES

26

Primero, como ξ es una raíz exacta de la ecuación f (x) = 0, se tiene que f (ξ) = 0 ⇐⇒ ξ = θ (ξ). Además, para cualquier k = 0, 1, 2, ... se tiene xk+1 = θ (xk ), desde aquí xk+1 − ξ = θ (xk ) − θ (ξ) Como θ es continua y diferenciable en I, por el teorema del valor medio, existe ck entre xk y ξ, tal que xk+1 − ξ = θ (xk ) − θ (ξ) = θ0 (ck ) (xk − ξ)

k = 0, 1, 2, ...

Luego, |xk+1 − ξ| = |θ0 (ck ) ||xk − ξ| ≤ M|xk − ξ| < |xk − ξ|

k = 0, 1, 2, ... (2.5)

Es decir, la distancia entre xk+1 y ξ es menor que la distancia entre xk y ξ, como I está centrado en ξ, vemos que si xk ∈ I entonces xk+1 ∈ I. Luego, si x0 ∈ I, claramente xk ∈ I para todo k = 1, 2, .... Para probar la segunda parte, l´ımk→∞ xk = ξ, desde (2.5) vemos que: |x1 − ξ| ≤ M|x0 − ξ| |x2 − ξ| ≤ M|x1 − ξ| ≤ M 2 |x0 − ξ| .. .

|xk − ξ| ≤ M k |x0 − ξ| de donde l´ımk→∞ |xk − ξ| = 0, pues 0 ≤ M < 1. Por lo tanto, l´ımk→∞ xk = ξ. Ejemplo 2.5 Sea la función f (x) = x2 + 2x − 10 cuya raíz es ξ ≈ 2,3166. Dadas las funciones iteración θ1 (x) = 5 − y θ2 (x) =

10 , x+2

x2 2 x 6= −2

Observe que |θ01 (x) | = | − x| = |x| < 1 ⇐⇒ x ∈ h−1, 1i Luego, no existe un intervalo I centrado en ξ tal que |θ01 (x) | < 1 para todo x ∈ I. El teorema 2.2 no afirma nada con respecto de la convergencia de la

CAPÍTULO 2. CEROS REALES DE FUNCIONES REALES

27

sucesión {xk } generada por xk+1 = θ1 (xk ), pues θ1 no cumple la hipótesis, el método del punto fijo puede convergir o no cuando se utilice θ1 como función iteración. Por otro lado, si usamos la función iteración θ2 , la situación es diferente. Observe que ¯ ¯ D E D√ E √ ¯ 10 ¯ 0 ¯ < 1 ⇐⇒ x ∈ −∞; − 10 − 2 ∪ ¯ |θ2 (x)| = ¯ 10 − 2; +∞ ≈ (x + 2)2 ¯ h−∞; −5,1622i ∪ h1,1622; +∞i Luego, existe un intervalo I centrado en ξ, tal que |θ02 (x)| < 1 para todo x ∈ I. En este caso, el teorema asegura la convergencia del MPF tomando x0 ∈ I. Ejercicio 2.4 Analice el caso para la función f (x) = x2 + 2x − 10, cuya raíz es ξ ≈ 2,3166, cuando se usa como función iteración: θ3 (x) =

10 − 2, x 6= 0 x

Además, analice el otro caso, cuando consideramos como raíz a ξ ≈ −5,1622. Algoritmo 2.2 (Punto Fijo) Considere la ecuación f (x) = 0 y la ecuación equivalente x = θ (x). Supongamos que θ ya es conocida explícitamente y las hipótesis suficientes del teorema 2.2 son satisfechas. Dada una precisión deseada ε > 0, el algoritmo se detendrá cuando |f (xk )| < ε. La aproximación al punto fijo ξ (la raíz buscada) será x¯. 1. Datos iniciales: a) x0 , la aproximación inicial b) ε1 y ε2 las precisiones deseadas 2. Si |f (x0 )| < ε1 , hacer x¯ = x0 , finalizar el algoritmo 3. k = 1 4. x1 = θ (x0 ) 5. Si |f (x1 )| < ε1 o si |x1 − x0 | < ε2 , hacer x¯ = x1 , finalizar el algoritmo 6. x0 = x1 7. k = k + 1. Volver al paso 4.

CAPÍTULO 2. CEROS REALES DE FUNCIONES REALES

28

Laboratorio 2.4 En algún lenguaje de programación de su preferencia, implemente el MPF y ejecute el programa sobre el problema del ejemplo 2.5. Observe que el método diverge cuando usamos θ1 , pruebe con varios puntos iniciales e interprete los resultados. Tasa de Convergencia del Método de Punto Fijo Cuando vimos el método de la bisección, notamos que era posible obtener una cota inferior para el número de iteraciones a ser realizadas por el algoritmo. Pero eso no siempre es posible, así, necesitamos de algunos parámetros que nos indiquen con qué rapidez la sucesión generada por un algoritmo converge a la solución deseada. Eso nos permitirá calificar a un algoritmo como lento o rápido. © ª Definición 2.2 (Tasa o rapidez de convergencia) Sea la sucesión r(k) generada por algún algoritmo, de modo que l´ımk→∞ r(k) → r¯. Asumamos que r(k) 6= r¯ para todo k = 0, 1, 2, ..., la tasa de convergencia del algoritmo es el supremo P de los números no negativos p satisfaciendo: ° ° (k+1) °r − r¯° l´ım p = β < ∞ k→∞ kr (k) − r ¯k Si P = 1 y el radio de convergencia β < 1, la sucesión se dice que tiene una tasa de convergencia lineal (por lo menos lineal). Si P ≥ 1 y β = 0, la sucesión tiene tasa de convergencia súper lineal. Si P = 2 y β < ∞, entonces diremos que la sucesión tiene una tasa de convergencia cuadrática. Una natural de ver esta situación es la siguiente: supongamos ª © manera que r(k) es una sucesión generada por un algoritmo la cual converge a la solución r¯, donde el algoritmo presenta una tasa de convergencia lineal, entonces: ° (k+1) ° ° ° °r − r¯° ≈ β °r(k) − r¯°, ∀k ≥ k0 | | {z } {z } ek+1

ek

nos dice que el error ek+1 cometido en la iteración k+1 es menor (linealmente) que el error ek en la iteración k, cuando k0 es grande.

En el caso de una tasa de convergencia cuadrática, es claro que el error cometido en la iteración k + 1 es aproximadamente el cuadrado del error cometido en la iteración anterior. Esto indica que para valores grandes de k0 y en las proximidades de r¯, el error ek+1 disminuye considerablemente con respecto a ek : ° ° °2 ° (k+1) °r − r¯° ≈ β °r(k) − r¯° , ∀k ≥ k0 | {z } | {z } ek+1

e2k

CAPÍTULO 2. CEROS REALES DE FUNCIONES REALES

29

Así, un algoritmo con tasa de convergencia cuadrática convergirá con mayor rapidez hacia un punto de acumulación, que uno que posee tasa de convergencia súper lineal. Por otro lado, un algoritmo con tasa de convergencia súper lineal será más rápido que uno con tasa de convergencia lineal. Más adelante veremos que el Método Secante tiene una tasa de convergencia P = 1,618... Proposición 2.1 (Tasa de convergencia de MPF) Asumamos que la sucesión generada por el Método del Punto Fijo converge a ξ, la raíz de f , entonces la tasa de convergencia es por lo menos lineal. Prueba. Desde (2.5), tenemos |xk+1 − ξ| = |θ0 (ck ) ||xk − ξ|,

ck entre xk y ξ,

k = 0, 1, 2, ...

Luego, |xk+1 − ξ| = |θ0 (ck ) |, |xk − ξ|

ck entre xk y ξ,

k = 0, 1, 2, ...

Tomando límites, por la continuidad de θ y θ0 , tenemos ¯ ³ ´¯ |xk+1 − ξ| ¯ ¯ 0 l´ım = l´ım |θ (ck ) | = ¯θ0 l´ım ck ¯ = |θ0 (ξ)| = β ≤ M < 1 k→∞ |xk − ξ| k→∞ k→∞

Así, vemos que el método del punto fijo posee una tasa de convergencia lineal, por este motivo muchas veces se le considera lento. Sin embargo, el método de bisección también posee esa misma velocidad, tal como lo afirma el siguiente ejercicio. Ejercicio 2.5 Muestre que la tasa de convergencia del método de la bisección es lineal.

2.2.4.

El Método de Newton-Raphson

Cuando estudiamos el método del punto fijo, notamos que: 1. Una de las condiciones para la convergencia es que |θ0 (x) | ≤ M < 1, para todo x ∈ I, donde I es un intervalo centrado en la raíz. 2. La convergencia del método será más rápida cuanto menor sea |θ0 (ξ) |.

CAPÍTULO 2. CEROS REALES DE FUNCIONES REALES

30

La segunda afirmación se debe al siguiente hecho: cuando analizamos la tasa de convergencia del MPF, vimos que l´ım

|xk+1 − ξ| = |θ0 (ξ) | < 1 |xk − ξ|

Entonces, acelerar la convergencia del MPF se conseguiría escogiendo una función iteración de modo que θ0 (ξ) = 0, pues en este caso la tasa de convergencia sería por lo menos súper-lineal. Hacia este objetivo, dada la ecuación f (x) = 0 cuya raíz es ξ. Consideramos la forma general de la función iteración θ (x) = x + A (x) f (x), donde A (ξ) 6= 0, con la nueva condición θ0 (ξ) = 0. Así, θ (x) = x + A (x) f (x) =⇒ θ0 (x) = 1 + A0 (x) f (x) + A (x) f 0 (x)

Luego, θ0 (ξ) = 1 + A0 (ξ) f (ξ) + A (ξ) f 0 (ξ) =⇒ θ0 (ξ) = 1 + A (ξ) f 0 (ξ) de donde θ0 (ξ) = 0 ⇐⇒ 1 + A (ξ) f 0 (ξ) = 0 ⇐⇒ A (ξ) = −

1 , f 0 (ξ)

f 0 (ξ) 6= 0

Esto nos motiva a definir

1 (ξ) Entonces, dada la ecuación f (x) = 0, la función iteración deseada será de la forma: f (x) θ (x) = x − 0 (2.6) , f 0 (x) 6= 0 f (x) Observe que A (ξ) = −

f0

(f 0 (x))2 − f (x) f 00 (x) f (x) f 00 (x) = , θ (x) = 1 − (f 0 (x))2 (f 0 (x))2 0

f 0 (x) 6= 0

y como f (ξ) = 0, se tiene que θ0 (ξ) = 0. Por lo tanto, usando la función iteración definida en (2.6) obtenemos un caso particular del método del punto fijo, a éste se le denomina método de Newton. En resumen, dado x0 ∈ R, el Método de Newton consiste en construir una sucesión definida por la regla xk+1 = xk −

f (xk ) , f 0 (xk )

f 0 (xk ) 6= 0,

k = 0, 1, 2, ...

(2.7)

CAPÍTULO 2. CEROS REALES DE FUNCIONES REALES

31

Otro Enfoque del Método de Newton-Raphson Desde el punto de vista geométrico, el método de Newton puede ser visto como la solución de un problema difícil, mediante la sucesiva resolución de problemas fáciles. Es decir, dada una aproximación inicial xk ∈ R a la raíz buscada, el problema difícil será hallar una raíz de la ecuación no lineal f (x) = 0, mientras que el problema fácil asociado será resolver la ecuación Lk (x) = 0, donde L es una función lineal afín que es parecida, al menos localmente, a la función no lineal f en torno al punto xk . Así, sea el problema (difícil) que consiste en hallar una raíz de f (x) = 0 y x0 ∈ R una aproximación inicial. Por el teorema de Taylor, existe δ > 0 tal que f (x) ≈ L0 (x) = f (x0 ) + f 0 (x0 ) (x − x0 ) para todo x ∈ hx0 − δ, x0 + δi. Luego, denotando por x1 la solución de la ecuación lineal L0 (x) = 0 (2.8) y asumiendo que f 0 (x0 ) 6= 0, entonces L0 (x) = 0 si, y sólo si, f (x0 ) + f 0 (x0 ) (x − x0 ) = 0 de donde x1 = x0 −

f (x0 ) f 0 (x0 )

Esperamos que x1 sea una mejor aproximación que x0 a la solución de f (x) = 0. Este procedimiento puede ser repetido iterativamente, creándose una sucesión {xk }∞ k=0 , donde xk+1 = xk −

f (xk ) , f 0 (xk )

f 0 (xk ) 6= 0,

k = 0, 1, 2, ...

Bajo algunas hipótesis se consigue que l´ımk→∞ xk = ξ. Convergencia del Método de Newton-Raphson A continuación damos las condiciones bajo las cuales se asegura la convergencia del método de Newton. Teorema 2.3 Sean f , f 0 y f 00 continuas en un intervalo abierto I que contiene en su interior la raíz ξ de f (x) = 0, donde f 0 (ξ) 6= 0, entonces existe

CAPÍTULO 2. CEROS REALES DE FUNCIONES REALES

32

Figura 2.6: Iteraciones Newton ¯ la sucesión un intervalo abierto I¯ ⊂ I conteniendo la raíz ξ, tal que si x0 ∈ I, {xk } generada por xk+1 = xk −

f (xk ) , f 0 (xk )

f 0 (xk ) 6= 0,

k = 0, 1, 2, ...

convergirá para ξ. Prueba. Observe que el método de newton es en realidad un MPF con . Así, para probar la convergencia debemos función iteración θ (x) = x − ff0(x) (x) ¯ probar que existe un I ⊂ I, centrado en ξ, tal que: 1. θ y θ0 son continuas en I¯ 2. |θ0 (x) | ≤ M < 1, para todo x ∈ I¯ 00

(x)f (x) 0 , θ0 (x) = f(f Vemos que θ (x) = x− ff0(x) 0 (x))2 y por hipótesis f (ξ) 6= 0. Como (x) f 0 es continua en I, es posible obtener un intervalo abierto I1 ⊂ I, ξ ∈ I1 , tal que f 0 (x) 6= 0 para todo x ∈ I1 .

Así, en el intervalo I1 se tiene que f, f 0 y f 00 son continuas y f 0 (x) 6= 0, ∀x ∈ I1 . Por lo tanto, θ y θ0 son continuas en I1 . 00

(x)f (x) 0 0 Como θ0 (x) = f(f 0 (x))2 , θ es continua en I1 y θ (ξ) = 0, entonces es posible obtener otro intervalo abierto I2 ⊂ I1 , centrado en ξ, tal que |θ0 (x) | < 1 para todo x ∈ I2 .

CAPÍTULO 2. CEROS REALES DE FUNCIONES REALES

33

Por lo tanto, hemos encontrado un intervalo abierto I¯ = I2 centrado en ξ ¯ Así, si x0 ∈ I, ¯ la sucesión donde θ y θ0 son continuas y |θ0 (x) | < 1, ∀x ∈ I. {xk } generada por la regla de correspondencia xk+1 = xk − f (xk ) /f 0 (xk ) converge hacia la raíz ξ de f (x) = 0. En pocas palabras, lo que dice el teorema anterior es que el método de Newton converge sólo si el punto inicial es tomado lo suficientemente próximo de la solución ξ, esta propiedad se conoce como convergencia local, la cual es una desventaja. Más adelante veremos que, en cierta forma, ese defecto se compensa con la rapidez con que converge el método. Algoritmo 2.3 (Newton-Raphson) Dada la ecuación f (x) = 0. Suponga que las hipótesis de suficiencia dadas en el teorema 2.3 son satisfechas. El algoritmo otorgará una aproximación x¯ a las raíz ξ. Paso inicial: Sea x0 una aproximaxión inicial de ξ y ε1 , ε2 > 0 las precisiones deseadas. Paso principal: Si |f (x0 ) | < ε1 , hacer x¯ = x0 y finalizar. Caso contrario, hacer: 1. k = 1 2. x1 = x0 − f (x0 ) /f 0 (x0 )

3. Si |f (x1 ) | < ε1 o si |x1 − x0 | < ε2 , hacer x¯ = x1 y finalizar. Caso contrario, hacer a) x0 = x1 b) k = k + 1. Volver al paso 2. Laboratorio 2.5 En algún lenguaje de programación de su preferencia, implemente el algoritmo de Newton y experiméntelo con diversos ejemplares. Compare sus resultados con los métodos anteriormente estudiados. Tasa de Convergencia del Método de Newton-Raphson Dada ξ una raíz exacta de la ecuación f (x) = 0. Sea {xk } la sucesión generada por el método de newton, tal que l´ımk→∞ xk = ξ. Debido a que el método de Newton es un caso particular de MPF, entonces debe tener por lo menos una tasa de convergencia lineal. Nosotros mostraremos que es mucho más que eso, debido a la condición θ0 (ξ) = 0, el método alcanzará una tasa cuadrática, es decir, en la definición 2.2, P = 2.

CAPÍTULO 2. CEROS REALES DE FUNCIONES REALES

34

Supongamos además que se satisfacen todas las hipótesis del teorema 2.3. Observe que f (xk ) f (xk ) f (xk ) xk+1 = xk − 0 =⇒ xk+1 −ξ = xk −ξ − 0 =⇒ ek+1 = ek − 0 f (xk ) f (xk ) f (xk ) donde ek = xk − ξ. Usando la serie de Taylor para f en torno al punto xk : 1 f (x) = f (xk ) + f 0 (xk ) (x − xk ) + f 00 (ck ) (x − xk )2 , ck está entre x y ck 2 Para x = ξ, tenemos 1 0 = f (xk ) + f 0 (xk ) (ξ − xk ) + f 00 (ck ) (ξ − xk )2 2 de donde 1 f (xk ) = f 0 (xk ) (xk − ξ) − f 00 (ck ) (xk − ξ)2 2 0 Dividiendo entre f (xk ), obtenemos f (xk ) f 00 (ck ) (xk − ξ)2 = (xk − ξ) − f 0 (xk ) 2f 0 (xk ) = ek −

f 00 (ck ) e2k 2f 0 (xk )

Luego, f (xk ) f 00 (ck ) e2k = ek − 0 = ek+1 0 2f (xk ) f (xk ) de donde

f 00 (ck ) ek+1 = 0 e2k 2f (xk ) Después de unos cálculos, vemos que θ00 (x) =

(2.9)

¡ ¢0 (f 0 (x))2 (f 0 (x) f 00 (x) + f (x) f 000 (x)) − (f (x) f 00 (x)) (f 0 (x))2 (f 0 (x))4

de donde θ00 es continua en ξ y f 00 (ξ) θ (ξ) = 0 f (ξ) 00

Llevando (2.9) al límite, tenemos |ek+1 | = k→∞ |ek |2 l´ım

|f 00 (ck ) | k→∞ |2f 0 (xk ) | l´ım

1 |f 00 (l´ımk→∞ ck ) | 1 |f 00 (ξ) | 1 = = |θ00 (ξ) | = β 0 0 2 |f (l´ımk→∞ xk ) | 2 |f (ξ) | 2 Por lo tanto, el método de Newton tiene una tasa de convergencia cuadrática. =

CAPÍTULO 2. CEROS REALES DE FUNCIONES REALES

2.2.5.

35

Método Secante

Uno de los inconvenientes en el método de Newton es el cálculo de la derivada de f en cada iteración, es decir, evaluar f 0 (xk ). Una manera de enfrentar esto es considerar una aproximación de la derivada: f 0 (xk ) ≈

f (xk ) − f (xk−1 ) , xk − xk−1

xk ≈ xk−1

donde xk y xk−1 son dos estimativas de la raíz exacta ξ. Para este caso, la función iteración Newton queda establecida como θ (xk ) = xk −

f (xk ) f (xk )−f (xk−1 ) xk −xk−1

= xk −

f (xk ) (xk − xk−1 ) f (xk ) − f (xk−1 )

O si lo prefiere, dadas dos estimativas iniciales para la raíz ξ de la ecuación f (x) = 0, x0 y x1 , el método secante consiste en construir una sucesión {xk } definida por la regla xk+1 = xk −

f (xk ) (xk − xk−1 ) , f (xk ) − f (xk−1 )

k = 1, 2, 3, ...

(2.10)

Las condiciones para la convergencia del método secante son prácticamente las mismas que las del método de Newton. Si bien la dificultad del cálculo explícito de la derivada fue evitada, lamentablemente el precio que se paga por esto es la disminución en la tasa de convergencia, que es súper-lineal. Una interpretación gráfica puede ser vista en la figura 2.7.

Figura 2.7: Método secante

CAPÍTULO 2. CEROS REALES DE FUNCIONES REALES

36

Convergencia del Método Secante Al igual que el Método de Newton, la convergencia del Método Secante está asegurada cuando los puntos iniciales x0 y x1 fueron tomados lo suficientemente próximos de la raíz buscada ξ. El teorema a seguir formaliza esta afirmación. Teorema 2.4 Sea ξ una raíz de f (x) = 0 y suponga que f ∈ C 2 en una vecindad en torno de ξ donde f 0 (ξ) 6= 0. Si los puntos iniciales x0 y x1 están suficientemente próximos de ξ, entonces l´ımk→∞ xk = ξ. Prueba. Desde la ecuación (2.10), tenemos ξ − xk+1 = ξ − xk +

f (xk ) (xk − xk−1 ) f (xk ) − f (xk−1 )  f (ξ)−f (x

 = − (ξ − xk−1 ) (ξ − xk )   µ

ξ−xk

k ) − f (xk )−f (xk−1 ) xk −xk−1

ξ−xk−1 f (xk )−f (xk−1 ) xk −xk−1

f 00 (γ k ) = − (ξ − xk−1 ) (ξ − xk ) 2f 0 (η k )

¶

   

(2.11)

donde γ k está entre ξ,xk y xk−1 , mientras que η k está entre xk y xk−1 . Como f ∈ C 2 , existe una vecindad I = [ξ − δ, ξ + δ], con δ > 0, tal que f 0 y f 00 son continuas y f 0 (x) 6= 0 para todo x ∈ I. Por tanto, existe M > 0 tal que ¯ 00 ¯ ¯ f (γ k ) ¯ ¯ ¯ M = m´ax ¯ 0 x∈I 2f (η k ) ¯

Elijamos los puntos iniciales x0 , x1 ∈ I, de modo que

M |ξ − x0 | < 1 y M |ξ − x1 | < 1 Definamos ahora t = m´ax {M |ξ − x0 | , M |ξ − x1 |}, claramente t < 1. Así, por (2.11) tenemos ¯ ¶¯ µ 00 ¯ f (γ k ) ¯¯ ¯ M |ξ − x2 | = M ¯(ξ − x1 ) (ξ − x0 ) 2f 0 (η k ) ¯ = M 2 |ξ − x1 | |ξ − x0 | ≤ t2 ≤ t < 1 de donde |ξ − x2 | <

t = m´ax {|ξ − x0 | , |ξ − x1 |} < δ M

CAPÍTULO 2. CEROS REALES DE FUNCIONES REALES

37

Así, x2 ∈ I. Usando este argumento, tenemos por inducción que xk ∈ I, para k ≥ 2. En general se tiene M |ξ − x2 | = M 2 |ξ − x1 | |ξ − x0 | ≤ t2 M |ξ − x3 | = M 2 |ξ − x2 | |ξ − x1 | ≤ t3 M |ξ − x4 | = M 2 |ξ − x3 | |ξ − x2 | ≤ t5 .. . M |ξ − xk+1 | = M 2 |ξ − xk | |ξ − xk−1 | ≤ tFk de donde Fk es el k-ésimo término de la sucesión de Fibonacci. Por lo tanto, l´ımk→∞ M |ξ − xk | = 0. Tasa de Convergencia del Método Secante A continuación analizaremos la tasa o rapidez de convergencia que el Método Secante posee. Debemos recalcar que algunos autores denominan a esta propiedad como el orden de convergencia. Teorema 2.5 Sea ξ una raíz de f (x) = 0. Asumamos que la sucesión generada por el Método Secante converge hacia ξ. Si f ∈ C 2 y √f 0 (ξ) 6= 0, entonces el Método secante tiene una tasa de convergencia p = 1+2 5 = 0,618... y ¶p−1 µ 00 ξ − xk+1 f (γ k ) l´ım p = − k→∞ (ξ − xk ) f 0 (η k ) Prueba. Desde (2.11) tenemos ξ − xk+1

µ

f 00 (γ k ) = − (ξ − xk−1 ) (ξ − xk ) 2f 0 (η k )

¶

(2.12)

donde γ k está entre ξ,xk y xk−1 , mientras que η k está entre xk y xk−1 . Resolviendo (2.12) obtenemos (ξ − xk ) =

1 pn rn A B K

(2.13)

´ 1 ³ √ √ f 00 (γ k ) ((ξ−x0 )K)r r−p 1+ 5 1− 5 donde K = − 2f , r = , A = y B = 0 (η ) , p = 2 2 (ξ−x1 )K k 1 ³ ´ r−p (ξ−x1 )K , la solución de (2.12) dada por (2.13) se puede comprobar ((ξ−x0 )K)p por sustitución directa. Por lo tanto, ξ − xk+1 p−1 rn B (r − p) p = K (ξ − xk )

CAPÍTULO 2. CEROS REALES DE FUNCIONES REALES y

38

µ 00 ¶p−1 ξ − xk+1 f (γ k ) p−1 rn l´ım ım K l´ım B (r − p) = − p = l´ k→∞ (ξ − xk ) k→∞ k→∞ 2f 0 (η k ) √

ξ−xk+1 1+ 5 Como se puede observar, l´ımk→∞ (ξ−x ≈ 1,61803398, p existe para p = 2 k) el cual en cierta forma nos recuerda al número áureo, 0,61803398..., concluimos entonces que la tasa o rapidez (orden) de convergencia del Método Secante es aproximadamente 1,618.

Como se puede observar, el Método Secante posee una tasa de convergencia un tanto inferior al Método de Newton, que tiene una tasa cuadrática. No obstante, la tasa de convergencia de 1,618 es superior a una lineal, y en la práctica esa disminución en la velocidad de convergencia se compensa con el hecho de no requerir el cálculo explícito de la derivada y su evaluación en cada iteración.. Laboratorio 2.6 En algún lenguaje de programación de su preferencia, implemente el método secante y experiméntelo con diversos ejemplares. Compare en la práctica sus resultados con los métodos anteriormente estudiados. Ejercicio 2.6 Experimente y compare los métodos estudiados en este capítulo, hallando una raíz de la ecuación f (x) = 0, donde f (x) = x − x ln (x) Ejercicio 2.7 Experimente y compare los métodos estudiados en este capítulo, hallando una raíz de la ecuación f (x) = 0, donde f (x) = ex − 4x2 Ejercicio 2.8 Sea x2 + x (ln (x) − 1) 2 Halle sus puntos críticos (puntos donde f 0 (x) = 0) usando un método iterativo estudiado en este capítulo. f (x) =

Ejercicio 2.9 Halle un punto donde la función f (x) = x2 + x3 ∗ (log(x) − 3) + 850 alcanza un mínimo sobre el intervalo [10, 20].

CAPÍTULO 2. CEROS REALES DE FUNCIONES REALES

39

Ejercicio 2.10 (ciclaje en el método de newton) Las iteraciones del Método de Newton entrarán en un ciclo ilimitado si xn+1 − a = − (xn − a) Esto a su vez sucede si f satisface x−a−

f (x) = − (x − a) f 0 (x)

La anterior expresión es una ecuación diferencial ordinaria separable de la forma f 0 (x) 1 =− f (x) 2(x − a) cuya solución es

p f (x) = sign(x − a) |x − a|

donde la función sign es la función signo, el cero de f es claramento x = a. Grafique f para el caso a = 2. Seguidamente, use un punto inicial x0 , donde x0 6= a, y ejecute el algoritmo de Newton. Ejercicio 2.11 Se considera la función F (x) = x5 + 2x. Mediante el Método de Newton, hallar el menor núumero positivo x (con tres decimales) para el cual F (x) = 4. Ejercicio 2.12 Probar, mediante el método de Newton, que la ecuación xn+1 = xn (2 − axn ) se puede utilizar para aproximar 1/a si x0 es una estimación inicial del recíproco de a. Nótese que este método de aproximar recíprocos utiliza sólo operaciones de suma y multiplicación. (Considerar f (x) = x1 − a). Pruebe para los casos: 1.

1 3

2.

1 11

Ejercicio 2.13 Pruebe que, al ser aplicado el método de Newton en la resolución de la ecuación ax+b = 0, donde a 6= 0, éste requiere sólo una iteración, sin importar qué punto inicial fue tomado. Ejercicio 2.14 En los casos siguientes, aplicar el método de Newton con la estimación inicial propuesta, explicar por qué falla el método.

CAPÍTULO 2. CEROS REALES DE FUNCIONES REALES

40

1. y = 2x3 − 6x2 + 6x − 1, donde x0 = 1 2. y = 4x3 − 12x2 + 12x − 3, donde x0 =

3 2

3. y = −x3 + 3x2 − x + 1, donde x0 = 1 √ 4. y = 3 x − 1, donde x0 = 2

2.3.

Comparación de los Métodos Iterativos

El esfuerzo computacional para la ejecución de cada uno de los métodos depende de varios factores, los más importantes son: 1. La complejidad de los cálculos, sobre todo para la derivada. 2. El número total de iteraciones 3. Condiciones para la convergencia El método de la bisección y el método de la posición falsa exigen pocas condiciones para garantizar la convergencia, el inconveniente está en que el número de iteraciones puede ser grande. Observe que su tasa de convergencia es lineal. Los métodos de punto fijo frecuentemente son más rápidos, pero a cambio exigen muchas hipótesis para la convergencia. El más rápido es el método de newton, pero requiere el cálculo de la derivada y demanda, al igual que los métodos de punto fijo, hipótesis rigurosas para su convergencia. El método secante puede ser práctico cuando el cálculo de la derivada es complicado, pero no es tan rápido como el método de Newton. Se puede concluir que la elección del método más eficiente depende de la ecuación que se intenta resolver. Cada método tiene sus ventajas y desventajas. Como un comentario adicional, después de llevar al computador cada uno de estos métodos y experimentarlos con diversos ejemplares, probablemente el estudiante halle que las diferencias de tiempo de ejecución, entre un programa y otro, sea insignificante cuando se aplica a la resolución de una ecuación, y ese afán por buscar el método más rápido parecería no tener sentido. Esa

CAPÍTULO 2. CEROS REALES DE FUNCIONES REALES

41

percepción es equivocada, pues estos métodos deben verse como subrutinas de otros métodos iterativos más sofisticados, para otro tipo de problemas, donde la pérdida de una fracción de segundo retrasaría el desempeño del método en su conjunto.

2.4.

Problemas

Problema 2.1 Considere las dos vigas de 30 m y 20 m cruzándose a una altura de 8 m del suelo, tal como se muestra en la figura 2.8. Determine el ancho del pasadizo, H.

Figura 2.8: Problema de las vigas

Problema 2.2 La concentración c de una bacteria contaminante en un lago decrece según la expresión c(t) = 80e−2t + 20e−0,5t , siendo t el tiempo en horas. Determinar el tiempo que se necesita para que el número de bacterias se reduzca a 7. Problema 2.3 Una determinada sustancia se desintegra según la ecuación A(t) = P e−0,0248t , donde P es la cantidad inicial en el tiempo t = 0 y A la cantidad resultante después de t años. Si inicialmente se depositan 500 miligramos de dicha sustancia, ¿cuánto tiempo debe transcurrir para que quede el 1 por ciento de ésta? Utilizar el Método de Newton. Problema 2.4 Una medicina administrada a un paciente produce una concentración en la sangre dada por c(t) = Ate−t/3 mg / ml, t horas después de que se hayan administrado A unidades. La máxima concentración sin peligro es de 1 mg / ml, y a esta cantidad se le denomina concentración de seguridad.

CAPÍTULO 2. CEROS REALES DE FUNCIONES REALES

42

1. ¿Qué cantidad debe ser inyectada para alcanzar como máximo esta concentración de seguridad? ¿Cuándo se alcanza este máximo? 2. Una cantidad adicional se debe administrar al paciente cuando la concentración baja a 0,25 mg / ml. Determínese con un error menor de 1 minuto cuándo debe ponerse esta segunda inyección. Problema 2.5 El crecimiento de poblaciones grandes puede modelarse en períodos cortos suponiendo que el crecimiento de la población es una función continua en t mediante una ecuación diferencial cuya solución es: N(t) = N0 eλt +

¢ v ¡ λt e −1 λ

donde N(t) es el número de individuos en el tiempo t (medido en años), λ es la razón de natalidad, N0 es la población inicial y v es un razón constante de inmigración, que se mide en número de inmigrantes al año. Supóngase que una población dada tiene un millón de individuos inicialmente y una inmigración de 400, 000 individuos al año. Se observa que al final del primer año la población es de 1506000 individuos. Se pide: 1. Determinar la tasa de natalidad. 2. Hacer una previsión de la población al cabo de tres años.

Capítulo 3 Resolución de Sistemas Lineales Sistemas de ecuaciones lineales aparecen en la resolución de diversos problemas de la vida real, frecuentemente, algunos de esos problemas involucran un gran número de variables y ecuaciones. Claramente, el método basado en la sustitución de variables, empleado para resolver pequeños sistemas de dos o tres variables, resulta obsoleto en estos casos, por lo que es necesario conocer métodos numéricos especializados para hacer frente a estas situaciones. Estos métodos numéricos no podrían trabajar sin la intervención de un computador, debido a la cantidad de cálculo envuelto en ese procedimiento. Todo esto sin lugar a dudas torna importante el estudio de sistemas de ecuaciones lineales desde el punto de vista numérico. Pero la resolución de sistemas de ecuaciones lineales no sólo tiene importancia propia, sino que constituye también una herramienta indispensable en la resolución de sistemas de ecuaciones no lineales. Más aún, es la base para la resolución de muchos otros problemas que surgen en diferentes áreas. Considere el siguiente sistema de ecuaciones lineales compuesto de dos ecuaciones y dos variables: ½ a1,1 x1 + a1,2 x2 = b1 a2,1 x1 + a2,2 x2 = b2

En sistemas simples como éste, podemos notar lo siguiente: 1. Solución única ½ 3x1 + 2x2 = 10 x1 − x2 = 5 2. Soluciones infinitas ½ 3x1 − 3x2 = 15 x1 − x2 = 5

=⇒

=⇒

43

· ¸ · ¸ 4 x¯1 = −1 x¯2

¸ · ¸ · 5+t x¯1 , = t x¯2

t∈R

CAPÍTULO 3. RESOLUCIÓN DE SISTEMAS LINEALES 3. Ninguna solución

½

44

2x1 − 2x2 = 8 x1 − x2 = 5

Gráficamente, cada ecuación representa a una recta, un punto que satisface ambas ecuaciones al mismo tiempo deberá estar en la intersección de ambas rectas (figura 3.1).

Figura 3.1: Representación gráfica de sistemas de ecuaciones lineales Por otro lado, cuando tenemos 3 variables, cada ecuación representaría un plano, un punto que satisface tres ecuaciones de un sistema simultáneamente estará en la intersección de tres planos. Para el caso n-dimensional, la situación es análoga, pero el conjunto definido por cada ecuación se denomina n-hiperplano.

3.1.

Aspectos Teóricos

Para un caso general, un sistema de m es la formulación matemática siguiente:   a1,1 x1 + a1,2 x2 + · · ·    a2,1 x1 + a2,2 x2 + · · · .. .. .. .. . .  . . . . .    a x + a x + ··· m,1 1

m,2 2

ecuaciones lineales y n variables + a1,n xn = b1 + a2,n xn = b2 .. .. .. .. . . . . + am,n xn = bm

(3.1)

Usando la notación matricial, el mismo sistema puede ser visto por Ax = b

(3.2)

CAPÍTULO 3. RESOLUCIÓN DE SISTEMAS LINEALES donde A ∈ Rm×n , definidas por:  a1,1  a2,1  A =  ..  . am,1

45

b ∈ Rm y x ∈ Rn es el vector de las variables o incógnitas, 

 a1,2 · · · a1,n a2,2 · · · a2,n   .. ..  , ... . .  am,2 · · · am,n

 b1  b2    b =  ..  . bm

y



 x1  x2    x =  ..  . xn

£ ¤t Una solución para el sistema (3.1) es un vector xˆ = xˆ1 · · · xˆn ∈ Rn el cual satisface simultáneamente cada una de las m ecuaciones que conforman el sistema. En un sistema lineal puede suceder lo siguiente: 1. El sistema lineal tiene única solución 2. El sistema lineal posee infinitas soluciones 3. El sistema lineal no posee soluciones Cuando n = m y A es inversible en (3.2), entonces la solución es única, en consecuencia, dicha solución puede ser obtenida haciendo x = A−1 b. No obstante, se debe advertir que este procedimiento tiene sólo un valor teórico, pues desde el punto de vista computacional es considerado costoso por el excesivo número de operaciones involucradas para calcular la inversa de A. En contraste, existen métodos más apropiados que no requieren el cálculo explícito de A−1 , como los que veremos más adelante.

3.1.1.

Conjunto Imagen de una Matriz A

Dada una matriz A ∈ Rm×n , definimos la imagen de A, denotado por Im (A), como el conjunto: Im (A) = {y ∈ Rm : ∃x ∈ Rn , y = Ax} Algunas veces a este conjunto se le donomina espacio imagen de A, debido a lo siguiente: Ejercicio 3.1 Muestre que Im (A) es un subespacio vectorial de Rm . Desde el punto de vista de las columnas de A ∈ Rm×n , note que Ax =

n X j=1

aj xj

CAPÍTULO 3. RESOLUCIÓN DE SISTEMAS LINEALES

46

donde aj representa la j-columna de A y xj es la j-componente del vector x ∈ Rn . Luego, resolver el sistema Ax = b significa obtener escalares x1 , ..., xn que permitan escribir b ∈ Rm como una combinación lineal de las columnas de A, es decir: b = a1 x1 + ... + an xn Por esta razón, al conjunto Im (A) también se le conoce como espacio columna de A. En consecuencia, la dimensión de la imagen de A, denotada por dim (Im (A)), está determinada por el máximo número de columnas linealmente independientes de A.

3.1.2.

Rango de una Matriz A

El rango de la matriz A ∈ Rm×n , denotado por rango (A), es definido por rango (A) = dim (Im (A)) Una manera práctica de encontrar el rango (A) consiste en determinar el mayor número de columnas linealmente independientes de A. En Álgebra Lineal se conoce que rango (A) = rango (At ), por consiguiente, para encontrar el rango (A) basta determinar también el mayor número de filas linealmente independientes de A, esto último suele ser de gran ayuda en muchos casos. Con respecto a la resolución del sistema de ecuaciones lineales Ax = b, podemos destacar lo siguiente: 1. Si m = n, puede pasar lo siguiente: a) Si rango (A) = n, el sistema Ax = b tiene solución única, pues A es inversible. En estas condiciones, se dice que el sistema es compatible y determinado. b) Si rango (A) < n, puede pasar lo siguiente 1) Si b ∈ Im (A), el sistema Ax = b admitirá infinitas soluciones. En estas condiciones, se dice que el sistema es compatible e indeterminado. 2) Si b 6∈ Im (A), el sistema no admite solución. En estas condiciones, se dice que el sistema es incompatible. 2. Si m > n, aunque rango (A) = n, el sistema puede no tener solución, debido a que es muy frecuente que b 6∈ Im (A).

CAPÍTULO 3. RESOLUCIÓN DE SISTEMAS LINEALES Matriz A ∈ Rm×n Rango completo Rango deficiente

mn rango (A) = n b ∈ Im (A) so lu c ió n ú n ica b 6∈ Im (A) in co m p a tib le

infinitas soluc. incompatible

Cuadro 3.1: Soluciones de un sistema de ecuaciones lineales 3. Si m < n, es muy probable que el sistema sea compatible e indeterminado si b ∈ Im (A). En el caso que b 6∈ Im (A), el sistema será incompatible. El cuadro 3.1 esquematiza las soluciones de un sistema de ecuaciones lineales.

3.2.

Métodos Numéricos para Resolver Sistemas de Ecuaciones Lineales

Los métodos numéricos que estudiaremos en este capítulo requieren que el sistema (3.1) esté constituido de n filas y n columnas, es decir, n = m. Al final de este capítulo, en el problema 3.4, será discutida una estrategia para enfrentar sistemas lineales donde m < n. Por otro lado, sistemas donde m > n se denominan sobredeterminados, ellos serán discutidos en el siguiente capítulo. Los métodos para resolver sistemas de ecuaciones lineales de n×n pueden ser de dos tipos, directos e iterativos. Directos: Si la solución existe, otorgan la solución exacta del sistema lineal después de un número finito de operaciones, excepto errores de redondeo. Iterativos: Dada una aproximación inicial x0 , generan una sucesión de vectores {xk }∞ k=0 . Si la solución existe, bajo ciertas condiciones, esta sucesión converge a la solución. A simple vista, parece ser que la elección de un método directo es la más adecuada. No obstante, como veremos más adelante, los métodos directos tienen un alto costo computacional, requiriendo alrededor de n3 operaciones elementales para resolver un sistema lineal de n×n. En contraste, los métodos iterativos parecen menos atractivos, sin embargo, sólo requieren alrededor de

CAPÍTULO 3. RESOLUCIÓN DE SISTEMAS LINEALES

48

n2 operaciones elementales por cada iteración1 , eso los torna más viables en la resolución de sistemas de ecuaciones lineales de grandes dimensiones2 . En las siguientes secciones analizaremos algunos de los métodos más importantes para la resolución numérica de sistemas de ecuaciones lineales.

3.3.

Método de Cramer

La “Regla de Cramer” puede ser enfocada como un método directo. Proviene de un teorema en álgebra lineal, mediante el cual se puede obtener la solución de un sistema lineal de ecuaciones en términos de determinantes. Recibe este nombre en honor a Gabriel Cramer (1704 - 1752). Si Ax = b es un sistema de ecuaciones lineales, donde A ∈ Rn×n es inversible y b ∈ Rn es un vector columna, entonces la solución del sistema se calcula así: det (Aj ) xj = j = 1, ..., n (3.3) det (A) donde Aj es la matriz que resulta de reemplazar la j-columna de A por b. Ejercicio 3.2 El determinante de una matriz real cuadrada de orden n está definido por n X (−1)j+1 a1,j det (A [1, j]) (3.4) det (A) = j=1

donde a1,j es el elemento en la 1-fila y j-columna de la matriz A, mientras que A [i, j] es la submatriz obtenida al eliminar la 1-fila y la j-columna de A. Muestre que el número de multiplicaciones necesarias para hallar det (A) es aproximadamente n! Es posible calcular el determinante de una matriz sin utilizar directamente la definición dada en (3.4), pero eso tomará por lo menos alrededor de n3 operaciones elementales. Por lo tanto, para fines prácticos, métodos como el de Cramer deben estar fuera de nuestro interés, debido a que resulta muy costoso en la resolución numérica de sistemas de ecuaciones lineales3 . 1

Pero pueden requerir muchas iteraciones. Frecuentemente, cuando n > 1000. 3 La regla de Cramer tiene un valor teórico, pues se utiliza en la demostración de muchas propiedades. 2

CAPÍTULO 3. RESOLUCIÓN DE SISTEMAS LINEALES

3.4.

49

Método de Gauss

Este método numérico es el más conocido y encuadra dentro de los métodos directos. Dado el sistema lineal de n × n   a1,1 x1 + a1,2 x2 +    a2,1 x1 + a2,2 x2 + .. .. .. ..  . . . .    a x + a x + n,1 1

n,2 2

· · · + a1,n xn = b1 · · · + a2,n xn = b2 .. .. .. . . . .. . . . . · · · + an,n xn = bn

(3.5)

El método de eliminación de Gauss consiste en transformar el sistema (3.5), de un modo equivalente, a un sistema triangular superior:  a01,1 x1 + a01,2 x2 + · · · + a01,n xn = b01     a02,2 x2 + · · · + a02,n xn = b02 (3.6) .. .. .. . . . ..  . . . .    a0n,n xn = b0n

El beneficio de esto es que se puede resolver el sistema triangular (3.6) de modo eficiente, es así que de la última ecuación de (3.6) tenemos xn =

b0n a0n,n

Luego, xn−1 puede ser obtenido mediante xn−1

b0n−1 − a0n−1,n xn = a0n−1,n−1

Y así sucesivamente, obtenemos xn−2 , xn−3 , ..., x2 , finalmente b01 − a01,2 x2 − a01,3 x3 − · · · − a01,n xn x1 = a01,1

3.4.1.

Resolviendo un Sistema Triangular

El siguiente algoritmo resuelve un sistema de ecuaciones lineales de orden n, el cual ya está en la forma triangular superior.

CAPÍTULO 3. RESOLUCIÓN DE SISTEMAS LINEALES

50

Algoritmo 3.1 (Resolución de un sistema triangular superior) Dado un sistema triangular superior Ax = b de orden n, con elementos de A sobre la diagonal no nulos. Los valores de las variables xn , xn−1 , ..., x2 , x1 son obtenidos mediante: xn = bn /an,n Para k = (n − 1) , ..., 1 s=0 Para j = (k + 1) , ..., n s = s + ak,j xj xk = (bk − s) /ak,k Ejercicio 3.3 Análogamente al algoritmo 3.1, diseñar un algoritmo que resuelva un sistema triangular inferior de orden n. Ejercicio 3.4 ¿Cuántas operaciones elementales (sumas, restas, multiplicaciones, divisiones y comparaciones) son necesarias para la ejecución del algoritmo 3.1 y el algoritmo planteado en el ejercicio 3.3? Laboratorio 3.1 En algún lenguaje de programación, implementar el algoritmo 3.1 y el algoritmo planteado en el ejercicio 3.3. La conversión de un sistema de orden n a un sistema equivalente, y triangular superior, es posible en virtud al siguiente teorema del álgebra lineal: Teorema 3.1 Sea Ax = b un sistema de ecuaciones lineales. Aplicando sobre las ecuaciones de este sistema una sucesión de operaciones elementales: 1. Cambiar dos ecuaciones 2. Multiplicar una ecuación por una constante no nula 3. Adicionar un múltiplo de una ecuación a otra ecuación ¯ = ¯b, el cual es equivalente4 al sistema oriobtenemos un nuevo sistema Ax ginal Ax = b. 4

Equivalente, en el sentido que tienen las mismas soluciones

CAPÍTULO 3. RESOLUCIÓN DE SISTEMAS LINEALES

3.4.2.

51

Conversión de un Sistema de Ecuaciones Lineales a uno Triangular Superior

El siguiente algoritmo usa el Teorema 3.1 y convierte un sistema de ecuaciones lineales Ax = b en un sistema triangular superior equivalente. Algoritmo 3.2 Dado el sistema lineal Ax = b de n × n. Suponga que el elemento ak,k 6= 0 al inicio de la etapa k: Para k = 1, ..., n − 1 Para i = k + 1, ..., n mi,k = ai,k /ak,k ai,k = 0 Para j = k + 1, ..., n ai,j = ai,j − mi,k ak,j bi = bi − mi,k bk Definición 3.1 El número ai,k mi,k = , ak,k

i = 2, ..., n,

k = 1, ..., n − 1

introducido en el algoritmo 3.2, lo denominaremos (i, k)-multiplicador de la matriz A. Además, ak,k se denomina el k-ésimo pivote. Ejemplo 3.1 Considere el siguiente sistema de ecuaciones lineales Ax = b, donde     1 2 3 0 2    2 2 4 10 4    y b = b(0) =  A = A(0) =   1 1 −2 1   6  1 0 0 2 0 Aplicando el Algoritmo 3.2, para m2,1 = 2, m3,1 = 1 y m4,1 = 1, de  1 2 3  0 −2 −2 A(1) =   0 −1 −5 0 −2 −3

Para k = 2, m3,2 = 0,5,  1  0 A(2) =   0 0

k = 1, los multiplicadores respectivos son donde:    0 2  0  10  (1)    y b =  4  1  2 −2

m4,2 = 1 y

 2 3 0 −2 −2 10   0 −4 −4  0 −1 −8

y

b(2)



 2  0   =  4  −2

CAPÍTULO 3. RESOLUCIÓN DE SISTEMAS LINEALES Finalmente, para k = 3  1  0 A(3) =   0 0

tenemos m4,3 = 0,25 y  2 3 0 −2 −2 10   y 0 −4 −4  0 0 −7

b(3)

52



 2  0   =  4  −3

Resta sólo aplicar el Algoritmo 3.1 para resolver el sistema triangular superior A(3) x = b(3) , de donde obtuvimos   −6/7  25/7   x=  −10/7  3/7

Un inconveniente que puede suceder en la aplicación de los algoritmos 3.2 y 3.1, es el cálculo del multiplicador mi,k , pues se necesita que ak,k 6= 0 en cada iteración. Pero el simple hecho que ak,k sea pequeño puede ocasionar que el multiplicador mi,k tome valores inmensamente grandes, lo que puede ocasionar a su vez el mal condicionamiento5 del sistema. Ejemplo 3.2 (Sistema mal condicionado) Utilice algún lenguaje de programación para resolver, usando los algoritmos 3.2 y 3.1, el siguiente sistema de ecuaciones: ½ 11x1 + 2x2 = 5 1016 x1 + 0,5x2 = 9

Después de aplicar consecutivamente los algoritmos, obtuvimos como solución: ¸ · 7,26691434300103 × 10−16 x= 2,5

Al verificar si realmente es solución, vemos que ¸ · ¸ · 5 5 6= b = Ax = 9 8,51691434300102

Este fenómeno se debe a que el pivote A1,1 = 11, cuando comparado con A2,1 = 1016 , es muy pequeño. Observe que en estas condiciones, el multiplicador a2,1 1016 m2,1 = = a1,1 11 es muy grande, lo que ocasionará impresiciones en los cálculos siguientes, pues el computador trabaja con precisión finita. 5

La matriz A tiende a ser no inversible en la práctica, aunque en teoría lo sea. Esto puede producir imprecisiones en los cálculos realizados por el computador.

CAPÍTULO 3. RESOLUCIÓN DE SISTEMAS LINEALES

53

¿Pero qué significa mal condicionamiento? Ese concepto es de gran ayuda en la resolución de sistemas lineales. A continuación formalizaremos ese término.

3.5.

Sensibilidad en Sistemas Lineales: Número de Condición

Considere el sistema Ax = b, donde A es inversible y b es un vector no nulo. Digamos que xˆ sea la única solución. Consideremos ahora una perturbación ∆b de b, el sistema lineal perturbado será (3.7)

Ay = b + ∆b

Observe que el sistema (3.7) también tiene única solución, digamos que ésta sea yˆ. Denotemos ∆ˆ x = yˆ − xˆ, en consecuencia, yˆ = xˆ + ∆ˆ x. Es natural esperar que cuando ∆b sea pequeño, entonces ∆ˆ x también lo sea. Para cuantificar el tamaño de vectores, usaremos la norma vectorial euclídea6 k·k. Así, la medida de ∆b relativa a b es k∆bk / kbk, mientras que la medida de ∆ˆ x relativa a xˆ es dada por k∆ˆ xk / kˆ xk. Por lo tanto, en términos más precisos, esperamos que cuando k∆bk / kbk sea pequeño, entonces k∆ˆ xk / kˆ xk también lo sea. Así, como yˆ es solución del sistema (3.7), entonces Aˆ y = b + ∆b =⇒ A (ˆ x + ∆ˆ x) = b +∆b =⇒ A (∆ˆ x) = ∆b =⇒ ∆ˆ x = A−1 (∆b) , Considerando la norma matricial inducida por la vectorial, kAk = m´axx6=0 kAxk kxk tenemos ° ° ° ° ∆ˆ x = A−1 (∆b) =⇒ k∆ˆ (3.8) xk = °A−1 (∆b)° ≤ °A−1 ° k∆bk

Por otro lado,

b = Aˆ x =⇒ kbk = kAˆ xk ≤ kAk kˆ xk =⇒ Por (3.8) y (3.9) tenemos

6

Si x ∈ Rn , kxk =

° k∆ˆ xk ° k∆bk ≤ °A−1 ° kAk kˆ xk kbk

p x21 + ... + x2n .

1 kAk ≤ kˆ xk kbk

(3.9)

CAPÍTULO 3. RESOLUCIÓN DE SISTEMAS LINEALES

54

Observe que, si k∆bk / kbk es pequeño y kA−1 k kAk es un número razonablemente pequeño, entonces k∆ˆ xk / kˆ xk debería ser también pequeño. Sin embargo, si kA−1 k kAk es extremadamente grande, a pesar que k∆bk / kbk sea pequeño, no hay garantía que k∆ˆ xk / kˆ xk también lo sea. En otras pa−1 labras, si kA k kAk es grande, podemos advertir que el sistema Ax = b puede ser susceptible a grandes alteraciones en la solución si b es ligeramente perturbado. Esto es crucial, pues el método de Gauss altera el vector b en cada iteración. Como se puede apreciar, el papel del número kA−1 k kAk juega un rol importante en la determinación de la estabilidad del sistema de ecuaciones lineales. A continuación lo formalizamos. Definición 3.2 (Número de condición) Dada la matriz A ∈ Rn×n inversible. El número de condición de A, denotado por cond (A), es ° ° cond (A) = °A−1 ° kAk

Si cond (A) es grande, diremos que la matriz es mal condicionada, caso contrario, ella será bien condicionada.

Ejercicio 3.5 Mostrar que, si k·k es una norma matricial inducida por la norma vectorial euclídea en Rn , cond (A) ≥ 1. Solución. Observe que la norma matricial usada es consistente, así: ° ° ° ° I = A−1 A =⇒ 1 = kIk = °A−1 A° ≤ °A−1 ° kAk = cond (A)

Observación 3.1 Algunas veces se utiliza también la siguiente expresión, rcond (A), la cual está definida por rcond (A) = 1/ cond (A) Claramente, si recond (A) es próximo a cero, la matriz A será mal condicionada.

3.6.

Estrategia con Pivotamiento Parcial

Como habíamos comentado anteriormente, uno de los inconvenientes en la aplicación del método de Gauss consistía en el cálculo de los multiplicadores, si el multiplicador mi,k = ai,k /ak,k es demasiado grande, éste podría ocasionar cálculos imprecisos en el sistema de cómputo. Algunas de las estrategias para evitar tal inconveniente son:

CAPÍTULO 3. RESOLUCIÓN DE SISTEMAS LINEALES

55

1. Estrategia con pivotamiento parcial 2. Estrategia con pivotamiento total Aquí veremos en qué consiste la estrategia de pivoteamiento parcial: En el inicio de la etapa k de la eliminación gaussiana, escoger como k−1 pivote el elemento de mayor módulo entre los coeficientes: ai,k , para k−1 i = k, k + 1, ..., n, donde ai,k representa la (i, k)-componente en la iteración k − 1. Permutar la fila k e i, si fuera necesario. Ejemplo 3.3 Considere un sistema de orden n, donde n = 4 y la iteración k = 2. Observe que   3 2 1 −1 5 £ (1) (1) ¤ 0 1 0 3 6  b A = 0 −3 −5 7 7 0 2 4 0 15 £ ¤ Note que A(1) b(1) representa la situación del sistema en la primera iteración. Ahora, para el inicio de la segunda iteración 1. Escoger el pivote: m´axj=2,3,4 |a1j,2 | = |a13,2 | = 3, entonces el pivote es −3 2. Cambiamos las líneas 2 y 3, de donde  3 2 £ (1) (1) ¤ 0 −3 b = A 0 1 0 2

obtenemos 1 −1 −5 7 0 3 4 0

 5 7  6 15

Observe que en este caso los multiplicadores respectivos serán: m3,2 = m4,2 =

1 1 =− −3 3 2 2 =− −3 3

Así, al escoger el mayor elemento en módulo entre los candidatos a pivote, se consigue que los multiplicadores m, en módulo, estén entre cero y uno, lo que evita la propagación de errores de redondeo.

CAPÍTULO 3. RESOLUCIÓN DE SISTEMAS LINEALES

56

La estrategia con pivoteamiento parcial no elimina del todo la acumulación de errores de redondeo, existe otra estrategia denominada estrategia con pivoteamiento completo. En contraste al pivoteamiento parcial, que busca el mejor pivote en una porción de cada columna en cada iteración, la estrategia de pivoteamiento completo analiza toda matriz. A pesar que en teoría esto elimina definitivamente las imprecisiones numéricas que puedan ocurrir, su uso no es cómun en la práctica pues requiere mucho esfuerzo computacional, es decir, requiere muchas operaciones elementales (comparaciones) para su ejecución. Problema 3.1 Investigue la estrategia de pivoteamiento completo.

3.7.

Descomposición LU

Esta estrategia consiste en descomponer la matriz A como el producto de dos matrices factores: A = LU, donde L es una matriz triangular inferior con unos sobre su diagonal, es decir   1 0 ··· 0  2,1 1 · · · 0    L =  .. .. ..  . .  . . .  . n,1 n,2 · · · 1

y U es una matriz triangular superior  u1,1 u1,2 · · · u1,n  0 u2,2 · · · u2,n  U =  .. .. .. ...  . . . 0 0 · · · un,n

    

De ese modo, resolver el sistema Ax = b, o sino LU x = b, es equivalente a resolver consecutivamente: 1. Ly = b 2. U x = y Una de las principales ventajas que tiene resolver sistemas mediante este método, en contraste con el método de Gauss, es la siguiente: supongamos que se quiere resolver los sistemas Ax = b y Ax = b0 , una vez obtenidos tales factores L y U , bastaría resolver consecutivamente los sistemas triangulares: Ly = b

Ux = y

(3.10)

CAPÍTULO 3. RESOLUCIÓN DE SISTEMAS LINEALES

57

y Ly = b0

(3.11)

Ux = y

Observe que los mismos factores L y U fueron usandos para resolver los sistemas (3.10) y (3.11), los cuales son fáciles de resolver usando el algoritmo 3.1 para sistemas triangulares superiores.

3.7.1.

Cálculo de los Factores L y U

El cálculo de las matrices L y U de modo que A = LU, está basado en los multiplicadores mi,j introducidos en la eliminación gaussiana (definión 3.1). Supongamos que tenemos la matriz   a01,1 a01,2 a01,3  0  0 0  a a a A0 =  2,1 2,2 2,3  =A a03,1 a03,2 a03,3

Primera iteración: Los respectivos multiplicadores están dados por m2,1 = a02,1 /a01,1 m3,1 = a03,1 /a01,1

para deshacernos de la variable x1 en las i-filas, para i = 2, 3, hacemos a11,j = a01,j

j = 1, 2, 3

a1i,j = a0i,j − mi,1 a01,j y obtenemos



i = 2, 3,

a11,1 a11,2 a11,3

 A1 =   0 0

Observe que A1 = M 0 A0 , donde 

j = 1, 2, 3 

 a12,2 a12,3   1 1 a3,2 a3,3

 1 0 0 M 0 =  −m2,1 1 0  −m3,1 0 1

CAPÍTULO 3. RESOLUCIÓN DE SISTEMAS LINEALES

58

Segunda iteración: Para eliminar x2 de las restantes líneas, calculamos el multiplicador m3,2 = a13,2 /a12,2 y hacemos a21,j = a11,j

j = 1, 2, 3

a22,j = a12,j

j = 2, 3

a23,j = = a13,j − m3,2 a12,j y obtenemos



a21,1 a21,2 a21,3

 A2 =   0 0

Observe que

j = 1, 2, 3 

 a22,2 a22,3   2 0 a3,3

A2 = M 1 A1 donde

 1 0 0 1 0  M1 =  0 0 −m3,2 1 

Además, A2 es triangular superior y podemos definir U = A2 . Luego: ¡ ¡ ¢ ¡ ¢ ¢ U = M 1 A1 = M 1 M 0 A0 = M 1 M 0 A0 = M 1 M 0 A

Es decir, y como

¡ ¢−1 ¡ 1 ¢−1 ¢−1 ¡ U = M0 U M A = M 1M 0

¡ 0 ¢−1 ¡ 1 ¢−1 M M

−1  −1 1 0 0 1 0 0 1 0  =  −m2,1 1 0   0 −m3,1 0 1 0 −m3,2 1      1 0 0 1 0 0 1 0 0 1 0  =  m2,1 1 0   0 1 0  =  m2,1 m3,1 0 1 m3,1 m3,2 1 0 m3,2 1 

concluimos que la matriz triangular inferior debería ser   1 0 0 1 0  L =  m2,1 m3,1 m3,2 1

de modo que A = LU .

CAPÍTULO 3. RESOLUCIÓN DE SISTEMAS LINEALES

59

Teorema 3.2 Dada una matriz cuadrada A de orden n, sea la matriz Ak constituida de las primeras k líneas y k columnas de A. Suponga que det (Ak ) 6= 0 para k = 1, ..., n − 1. Entonces, existe una única matriz triangular inferior L = [mi,j ], con mi,i = 1, para i = 1, ..., n, y una única matriz triangular superior U = [ui,j ] tales que A = LU. Además, det (A) = u1,1 u2,2 ...un,n . Ejercicio 3.6 Investigue la demostración del teorema 3.2. Ejemplo 3.4 Resolver el siguiente sistema de ecuaciones lineales utilizando la descomposición LU. x1 + 2x2 + 3x3 = 10 2x1 + 5x2 − x3 = 20 −x1 + 2x2 + x3 = 6 Observe que

Luego,



 1 2 3 A =  2 5 −1  −1 2 1 m2,1 = 2,

y

m3,2 = 4,

 10 b =  20  6 

m3,1 = −1,

y





 1 2 3 A1 =  0 1 −7  0 4 4

 1 2 3 A2 =  0 1 −7  0 0 32

Por lo tanto,        1 2 3 1 0 0 1 0 0 1 2 3 1 0  A2 =  2 1 0   0 1 −7  = LU A =  2 5 −1  =  m2,1 −1 2 1 −1 4 1 0 0 32 m3,1 m3,2 1

Puesto que Ax = b lo podemos expresar como (LU ) x = b, usando el algoritmo 3.1 y el ejercicio 3.3 para resolver sistemas triangulares, calculamos consecutivamente:   10 Ly = b =⇒ y =  0  16 y   3/2 Ux = y =⇒ x =  7/2  1/2

CAPÍTULO 3. RESOLUCIÓN DE SISTEMAS LINEALES

3.8.

60

Descomposición de Cholesky

Una estrategia que tiene suma importancia, sobre todo en optimización, es la descomposición de Cholesky. Definición 3.3 (Matriz definida positiva) Se dice que una matriz A ∈ Rn×n simétrica y de orden n es definida positiva, si xt Ax > 0, para todo x ∈ Rn y x 6= 0. Teorema 3.3 (Descomposición de Cholesky) Si A es una matriz de orden n, simétrica y definida positiva, entonces existe una única matriz triangular inferior de orden n y con diagonal positiva, tal que A = GGt Ejercicio 3.7 Investigue la demostración del teorema 3.3.

3.8.1.

Cálculo del factor de Cholesky

Dada una matriz simétrica y definida positiva:  a1,1 a1,2 · · · a1,n  a2,1 a2,2 · · · a2,n  A =  .. .. .. ...  . . . an,1 an,2 · · · an,n

El factor G será obtenido de la ecuación    g1,1 0 a1,1 a1,2 · · · a1,n  a2,1 a2,2 · · · a2,n   g2,1 g2,2     .. .. .. ..  =  .. ...  . . . .   . gn,1 gn,2 an,1 an,2 · · · an,n

matricial: ··· ··· ...

· · · gn,n

La columna 1:    g1,1 0 · · · a1,1 0  a2,1   g2,1 g2,2 · · · 0     ..  =  .. .. .. . ..  .   . . . an,1 gn,1 gn,2 · · · gn,n

luego

g1,1 =

√ a1,1

y

gj,1 =

0 0 .. .

    

aj,1 g1,1

g1,1 0 .. . 0

          

g1,1 g2,1 · · · gn,1 0 g2,2 · · · gn,2 .. .. .. ... . . . 0 0 · · · gn,n 

    =  

g1,1 g1,1 g2,1 g1,1 .. . gn,1 g1,1

para j = 2, ..., n

    

    

CAPÍTULO 3. RESOLUCIÓN DE SISTEMAS LINEALES La columna 2:    a1,2 g1,1 0 · · · 0  a2,2     g2,1 g2,2 · · · 0  a3,2     =  .. .. .. . ..  ..   . . .  .  gn,1 gn,2 · · · gn,n an,2

de donde





     

g1,2 g2,2 0 .. . 0





g1,1 g1,2 2 2 g2,1 + g2,2 g3,1 g2,1 + g3,2 g2,2 .. .

      =    

gn,1 g2,1 + gn,2 g2,2

g1,2 g1,1 = a1,2 =⇒ g1,2 = g2,1 = 2 2 g2,1 + g2,2 = a2,2 =⇒ g2,2 =

Como gj,1 g2,1 + gj,2 g2,2 = aj,2

a1,2 g1,1

q 2 a2,2 − g2,1

j = 3, ..., n

y los gj,1 , j = 1, ..., n, ya fueron calculados, entonces gj,2 =

(aj,2 − gj,1 g2,1 ) g2,2

j = 3, ..., n

Columna k: Los elementos de la columna k de G son: ¤t £ 0 · · · 0 gk,k gk+1,k · · · gn,k k = 2, ..., n

haciendo

tenemos de donde



ak,1 ak,2 .. .

      ak,k   ak+1,k  .  .. an,k



   g1,1 0 · · · 0    g2,1 g2,2 · · · 0    =  .. .. .. ...   . . .   gn,1 gn,2 · · · gn,n 

gk,1 gk,2 .. .



            gk,k     0   .   ..  0

2 2 2 ak,k = gk,1 + gk,2 + · · · + gk,k

gk,k = Y como



Ã

ak,k −

k−1 X i=1

2 gk,j

!1/2

aj,k = gj,1 gk,1 + gj,2 gk,2 + · · · + gj,k gk,k

j = k + 1, ..., n

61

      

CAPÍTULO 3. RESOLUCIÓN DE SISTEMAS LINEALES

62

y como los elementos gi,k , i = 1, ..., k − 1, ya fueron calculados, tenemos ³ ´ Pk−1 aj,k − i=1 gj,i gk,i j = k + 1, ..., n gj,k = gk,k Algoritmo 3.3 Sea A una matriz simétrica y definida positiva. Para k = 1, ..., n suma = 0; Para j = 1, ... (k − 1) 2 suma = suma + gk,j r = ak,k√− suma ..............(*) gk,k = r Para i = k + 1, ..., n suma = 0 Para j = 1, ..., k − 1 suma = suma + gi,j gk,j gi,k = (ai,k − suma) /gk,k Por lo general, ver si una matriz simétrica es definida positiva usando la definición es una tarea prácticamente imposible. Sin embargo, podemos usar el algoritmo de Cholesky para verificar si A es definida positiva. Si en (*) se tiene que r ≤ 0, entonces la descomposición no es posible y A no es definida positiva. Caso contrario, el algoritmo otorga la matriz triangular inferior G tal que A = GGt . La descomposición de Cholesky requiere alrededor de n3 /6 operaciones de multiplicación para la descomposición. Este número es aproximadamente la mitad del número de operaciones necesarias para la realización de la eliminación en la descomposición LU, pues requería n3 /3. Al igual que la descomposición LU, una vez conocido A = GGt , es posible resolver el sistema lineal asociado: Ax = b, del siguiente modo: como ¢ ¡ Ax = b ⇐⇒ GGt x = 0 entonces

1. Resolver Gy = b 2. Resolver Gt x = y

CAPÍTULO 3. RESOLUCIÓN DE SISTEMAS LINEALES

63

Laboratorio 3.2 En algún lenguaje de programación de su preferencia, implemente el algoritmo de Cholesky, de modo que: cuando la matriz A sea definida positiva, devuelva la matriz triangular inferior G. Caso contrario, debería emitir un mensaje anunciando que A no es definida positiva. Ejercicio 3.8 Usando el algoritmo de Cholesky, verifique que la matriz   27 12 10 24  12 9 16 29     10 16 42 70  24 29 70 137 es realmente definida positiva.

Ejercicio 3.9 Una duda común es, ¿podría ser definida positiva una matriz la cual tiene algunas componentes negativas? Verifique si la matriz   27 12 10 24  12 9 16 11     10 16 42 −2  24 11 −2 137

es o no definida positiva (use el algoritmo).

Ejercicio 3.10 Una matriz simétrica A ∈ Rn×n se llama semidefinida positiva, si xt Ax ≥ 0, para todo x ∈ Rn . Probar que, para cualquier matriz B ∈ Rm×n , la matriz C = B t B es semidefinida positiva. Solución: Observe que C ∈ Rm×m es simétrica, pues ¡ ¢t ¡ ¢t C t = BtB = Bt Bt = BtB = C

y

¢ ¡ xt Cx = xt B t B x = (Bx)t (Bx) = kBxk2 ≥ 0

donde k·k es la norma euclídea.

Ejercicio 3.11 Probar que, si A ∈ Rn×n es definida positiva, entonces A es inversible. Lo recíproco no siempre es válido, de un contraejemplo. Solución: Recordemos el teorema de álgebra lineal que nos decía: “Una matriz A ∈ Rn×n es inversible si, y sólo si, el sistema Ax = 0 tiene como única solución la trivial”. Ahora, con respecto al ejercicio, realizaremos una prueba por contradicción. Asumamos que A es definida positiva y no inversible. Entonces, el

CAPÍTULO 3. RESOLUCIÓN DE SISTEMAS LINEALES

64

sistema Ax = 0 tiene una solución no trivial, digamos que existe xˆ 6= 0 tal que Aˆ x = 0. Luego, 0 < xˆt Aˆ x = xˆt 0 = 0 Lo cual es una contradiccion. Esto completa la prueba. Para verificar la segunda parte del ejercicio, basta hacer notar que la matriz · ¸ −1 0 0 −1 es inversible pero no es definida positiva (en realidad es definida negativa).

3.9.

Métodos Iterativos

La idea de los métodos iterativos para resolver sistemas de ecuaciones está inspirada en el método de punto fijo. Dado el sistema lineal de ecuaciones Ax = b

(3.12)

donde A ∈ Rn×n y b ∈ Rn . Los métodos iterativos consisten en expresar el sistema (3.12) en una ecuación matricial equivalente: x = Dx + d

(3.13)

donde D ∈ Rn×n y d ∈ Rn . Note que en este caso θ (x) = Dx + d es una función iteración n-dimensional. Así, los métodos de este tipo construyen una sucesión {xk ∈ Rn }∞ k=0 definida por la regla xk+1 = θ (xk ) = Dxk + d Observe que, si l´ımk→∞ xk = x∗ , entonces x∗ es un punto fijo de θ. Bajo ciertas condiciones, x∗ debería ser también la solución del sistema original Ax = b. Método de Gauss-Jacobi Este método se caracteriza por transformar el sistema Ax = b a la forma x = Dx + d, del siguiente modo:

CAPÍTULO 3. RESOLUCIÓN DE SISTEMAS LINEALES

65

Sea el sistema original a1,1 x1 + a1,2 x2 + · · · + a1,n xn = b1 a2,1 x1 + a2,2 x2 + · · · + a2,n xn = b2 .. .. .. .. . . . .. .. .. . .. . . . . . . . an,1 x1 + an,2 x2 + · · · + an,n xn = bn Suponiendo que ai,i 6= 0 para i = 1, ..., n, despejamos las variables x1 , ..., xn de las n ecuaciones, respectivamente: x1 =

1 (b1 − a1,2 x2 − a1,3 x3 − ... − a1,n xn ) a1,1

1 (b2 − a2,1 x2 − a2,3 x3 − ... − a2,n xn ) a2,2 .. . 1 = (bn − an,1 x1 − an,2 x2 − ... − an,n−1 xn−1 ) an,n

x2 =

xn

De esta forma tenemos  b   1 0 − aa1,2 − aa1,3 · · · − aa1,n a1,1   1,1 1,1 1,1   x1  b2   a a a − a2,1 0 − a2,3 · · · − a2,n  x2   a2,2   2,2 2,2 2,2        ..  =  .  +  . .. .. .. ...  .   ..   .. . . .     bn  xn an,1 an,2 an,3 0 − an,n − an,n − an,n · · · an,n Es decir, la función iteración quedaría establecido por



       

x1 x2 .. . xn

    

x = θ (x) = Dx + d donde 

− aa1,2 1,1

− aa1,3 1,1

− aa1,n 1,1

0 ···  a2,1  − 0 − aa2,3 · · · − aa2,n  a2,2 2,2 2,2 D=  .. .. .. .. ...  . . . .  n,1 n,2 n,3 − aan,n − aan,n ··· 0 − aan,n



       

y

    d=   

b1 a1,1 b2 a2,2

.. . bn an,n

        

(3.14)

El método de Gauss-Jacobi consiste en lo siguiente: dado x0 , una aproximación inicial a la solución del sistema lineal Ax = b, obtener x1 , x2 , ..., xk , ...por medio de la relación recursiva xk+1 = Dxk + d

CAPÍTULO 3. RESOLUCIÓN DE SISTEMAS LINEALES

66

donde D y d están definidos en este caso por (3.14). O sino: ¢ 1 ¡ b1 − a1,2 xk2 − a1,3 xk3 − ... − a1,n xkn a1,1 ¢ 1 ¡ = b2 − a2,1 xk2 − a2,3 xk3 − ... − a2,n xkn a2,2 .. . ¢ 1 ¡ = bn − an,1 xk1 − an,2 xk2 − ... − an,n−1 xkn−1 an,n

xk+1 = 1 xk+1 2

xk+1 n

(3.15)

Una característica de los métodos iterativos, a diferencia de los métodos directos, es que sólo convergen si algunas hipótesis son satisfechas. Criterio de convergencia para el método de Gauss-Jacobi Dado el sistema lineal Ax = b, sea Pn j=1j6=i |ai,j | αi = y α = m´ax {αi } 1≤i≤n |ai,i | Si α