Métodos numéricos  

1. Introducción 2. Errores

o 2.1 Definiciones o 2.2 Dígitos significativos o 2.3 Propagación de errores o 2.4 Ejercicios adicionales

3. Aritmética de computadores o 3.1 Aritmética de punto fijo o 3.2 Números en punto flotante

3.2.1 Notación científica normalizada 3.2.2 Representación de los números en punto flotante

o 3.3 Aritmética de punto flotante 3.3.1 Números de máquina aproximados 3.3.2 Las operaciones básicas

o 3.4 Desbordamiento por exceso y desbordamiento por defecto o 3.5 Condicionamiento y estabilidad

4. Cálculo de raíces de ecuaciones o 4.1 Método de la bisección o 4.2 Método de las aproximaciones sucesivas o 4.3 Método de Newton o 4.4 Método de la secante o 4.5 Método de Steffensen o 4.6 Método de la falsa posición

5. Puntos fijos e iteración funcional 6. Resolución de sistemas de ecuaciones lineales

o 6.1 Métodos de resolución exacta 6.1.1 Sistemas   fáciles   de resolver 6.1.2 La factorización LU 6.1.3 Eliminación gaussiana básica 6.1.4 Método de Gauss-Jordan 6.1.5 Pivoteo

o 6.2 Métodos iterativos 6.2.1 Conceptos básicos 6.2.2 Método de Richardson 6.2.3 Método de Jacobi 6.2.4 Método de Gauss-Seidel

7. Interpolación o 7.1 Polinomios de interpolación de Lagrange o 7.2 Interpolación de splines o 7.3 Splines cúbicos

8. Integración numérica o 8.1 Integración vía interpolación polinomial o 8.2 Regla del trapecio o 8.3 Regla de Simpson

About this document ...

Wladimiro Diaz Villanueva

1. IntroducciónLa ciencia y la tecnología describen los fenómenos reales mediante modelos matemáticos. El estudio de estos modelos permite un conocimiento más profundo del fenómeno, así como de su evolución futura. La matemática aplicada es la rama de las matemáticas que se dedica a buscar y aplicar las herramientas más adecuadas a los problemas basados en estos modelos. Desafortunadamente, no siempre es posible aplicar métodos analíticos clásicos por diferentes razones:

No se adecúan al modelo concreto. Su aplicación resulta excesivamente compleja. La solución formal es tan complicada que hace imposible cualquier

interpretación posterior. Simplemente no existen métodos analíticos capaces de proporcionar

soluciones al problema.

En estos casos son útiles las técnicas numéricas, que mediante una labor de cálculo más o menos intensa, conducen a soluciones aproximadas que son siempre numérica. El importante esfuerzo de cálculo que implica la mayoría de estos métodos hace que su uso esté íntimamente ligado al empleo de computadores. De hecho, sin el desarrollo que se ha producido en el campo de la informática resultaría difícilmente imaginable el nivel actual de utilización de las técnicas numéricas en ámbitos cada día más diversos1.

2. ErroresEl concepto de error es consustancial con el cálculo numérico. En todos los problemas es fundamental hacer un seguimiento de los errores cometidos a fin de poder estimar el grado de aproximación de la solución que se obtiene.

Los errores asociados a todo cálculo numérico tienen su origen en dos grandes factores:

Aquellos que son inherentes a la formulación del problema. Los que son consecuencia del método empleado para encontrar la solución

del problema.

Dentro del grupo de los primeros, se incluyen aquellos en los que la definición matemática del problema es sólo una aproximación a la situación física real. Estos errores son normalmente despreciables; por ejemplo, el que se comete al obviar los efectos relativistas en la solución de un problema de mecánica clásica.

En aquellos casos en que estos errores no son realmente despreciables, nuestra solución será poco precisa independientemente de la precisión empleada para encontrar las soluciones numéricas.

Otra fuente de este tipo de errores tiene su origen en la imprecisión de los datos físicos: constantes físicas y datos empíricos. En el caso de errores en la medida de los datos empíricos y teniendo en cuenta su carácter generalmente aleatorio, su tratamiento analítico es especialmente complejo pero imprescindible para contrastar el resultado obtenido computacional-mente.

En lo que se refiere al segundo tipo de error (error computacional), tres son sus fuentes principales:

1.Equivocaciones en la realización de las operaciones (errores de bulto). Esta fuente de error es bien conocida por cualquiera que haya realizado cálculos manualmente o empleando una calculadora. El empleo de computadores ha reducido enormemente la probabilidad de que este tipo de errores se produzcan. Sin embargo, no es despreciable la probabilidad de que el programador cometa uno de estos errores (calculando correctamente el resultado erróneo). Más aún, la presencia de bugs no detectados en el compilador o en el software del sistema no es inusual. Cuando no resulta posible verificar que la solución calculada es razonablemente correcta, la probabilidad de que se haya cometido un error de bulto no puede ser ignorada. Sin embargo, no es esta la fuente de error que más nos va a preocupar.

2.El error causado por resolver el problema no como se ha formulado, sino mediante algún tipo de aproximación. Generalmente está causado por la sustitución de un infinito (sumatorio o integración) o un infinitesimal (diferenciación) por una aproximación finita. Algunos ejemplos son:

El cálculo de una función elemental (por ejemplo, Seno x) empleando sólo n términos de los infinitos que constituyen la expansión en serie de Taylor.

Aproximación de la integral de una función por una suma finita de los valores de la función, como la empleada en la regla del trapezoide.

Resolución de una ecuación diferencial reemplazando las derivadas por una aproximación (diferencias finitas).

Solución de la ecuación f(x) = 0 por el método de Newton-Raphson: proceso iterativo que, en general, converge sólo cuando el número de iteraciones tiende a infinito.

Denominaremos a este error, en todas sus formas, como error por truncamiento, ya que resulta de truncar un proceso infinito para obtener un proceso finito. Obviamente, estamos interesados en estimar, o al menos acotar, este error en cualquier procedimiento numérico.

3.Por último, la otra fuente de error de importancia es aquella que tiene su origen en el hecho de que los cálculos aritméticos no pueden realizarse con precisión ilimitada. Muchos números requieren infinitos decimales para ser representados correctamente, sin embargo, para operar con ellos es necesario redondearlos. Incluso en el caso en que un número pueda representarse exactamente, algunas operaciones aritméticas pueden dar lugar a la aparición de errores (las divisiones pueden producir números que deben ser redondeados y las multiplicaciones dar lugar a más dígitos de los que se pueden almacenar). El error que se introduce al redondear un número se denomina error de redondeo.

2.1 Definiciones

Ahora que disponemos de una idea correcta de qué es el error y de cual es su origen, podemos formalizar el concepto de error. Generalmente, no conocemos el valor de una cierta magnitud   y hemos de conformarnos con un valor aproximado x. Para estimar la magnitud de este error necesitamos dos definiciones básicas:

Error absoluto

de x: 

(1)

Error relativo

de x: 

  (2)

En la práctica, se emplea la expresión: 

  (3)

En general, no conocemos el valor de este error, ya que no es habitual disponer del valor exacto de la magnitud, sino sólo de una acotación de su valor, esto es, 

un número  , tal que: 

(4)

o bien: 

(5)

De acuerdo con este formalismo, tenemos que un numero se representará del siguiente modo: 

= (6)

= (7)

2.2 Dígitos significativosSea x un número real que, en general, tiene una representación decimal infinita. Podemos decir que x ha sido adecuadamente redondeado a un número con ddecimales, al que denominaremos x(d), si el error de redondeo,   es tal que: 

(8)

Ejemplo 1: Exprese el número x=35.47846 correctamente redondeado a cuatro (x(4)) y tres (x(3)) decimales. Calcular el error cometido.

Solución: en el primer caso obtenemos: 

x(4) = 35.4785  

=  

En el segundo caso, la aproximación correcta es: 

x(3) = 35.478  

=  

y no la siguiente: 

x(3) = 35.479  

=  

Es decir, no es correcto redondear por exceso cuando el dígito anterior es 5 y proviene de un acarreo previo. 

Otra forma de obtener el número de cifras significativas es mediante truncamiento, en donde simplemente se eliminan los dígitos de orden inferior. El error cometido en este caso es: 

(9)

y que, en general, conduce a peores resultados que el método anterior.

Ejemplo 2: Exprese el número x=35.47846 truncado a cuatro (x(4)) y tres (x(3)) decimales. Calcular el error cometido.

Solución: 

x(4) = 35.4784  

=  

x(3) = 35.478  

=  

2.3 Propagación de errores

Cuando se resuelve un problema matemático por métodos numéricos y aunque las operaciones se lleven a cabo exactamente, obtenemos una aproximación

numérica del resultado exacto. Es importante tratar de conocer el efecto que sobre el resultado final del problema tiene cada una de las operaciones realizadas.

Para estudiar como se propaga en error, veamos cual es el efecto que cada una de las operaciones básicas tiene sobre el error final cuando se aplican sobre dos

números   y   : 

= (10)

= (11)

= (12)

= (13)

Cuando el problema consiste en calcular el resultado y = f(x)tenemos la siguiente fórmula aproximada de propagación del error: 

(14)

En el caso más general, en que una función depende de más de una variable

(  ), la fórmula aproximada de propagación del error maximal es:

 (15)

Ejemplo 3: Determinar el error máximo cometido en el

cálculo y = x1 x22 para   y  .

Solución: El error cometido, de acuerdo con la ecuación (15), se puede calcular mediante: 

Sustituyendo valores, obtenemos: 

Por lo que el resultado final se debe expresar como:

Ejemplo 4: Sea el siguiente sistema de ecuaciones lineales: 

en donde  ; b = 1 / a y d = b - a¿Con qué exactitud podemos determinar el producto xy?

Solución: Primero resolveremos el sistema de ecuaciones por reducción: 

Ecuaciones que conducen a la siguiente expresión para el producto: 

 

(16)

Resolveremos ahora el problema por dos métodos. Primero, calcularemos el error asociado a cada una de las variables y los términos de la expresión anterior: 

Sustituyendo valores, obtenemos el siguiente resultado:

Una forma mucho más adecuada de resolver este problema consiste en sustituir en la expresión (16) los valores de b y d por sus correspondientes expresiones en función de a. Sustituyendo y operando, obtenemos que el producto y el error asociado vienen dados por: 

que, sustituyendo valores, conduce al resultado:

Si ambos resultados son correctos ¿Por qué el error es mucho menor en el segundo caso que en el primero? La respuesta es simple: en el segundo caso hemos eliminado operaciones intermedias, permitiendo que algunos errores se cancelen mutuamente. En general, cuanto menor sea el número de pasos intermedios que efectuemos para alcanzar la solución, menor será el error cometido.

2.4 Ejercicios adicionales1.

¿Con qué exactitud es necesario medir el radio de una esfera para que su volumen sea conocido con un error relativo menor de 0.01%? ¿Cuantos decimales es necesario emplear para el valor de  ?

Soluciones:  . El número   debe expresarse al menos con seis cifras decimales.

2.

Supongamos una barra de hierro de longitud l y sección 

rectangular   fija por uno de sus extremos. Si sobre el extremo libre aplicamos una fuerza Fperpendicular a la barra, la flexión s que ésta experimenta viene dada por la expresión: 

en donde E es una constante que depende sólo del material denominada módulo de Young. Conociendo que una fuerza de 140 Kp aplicada sobre 

una barra de 125 cm de longitud y sección cuadrada de 2.5 cm produce una flexión de 1.71 mm, calcular el módulo de Young y el intervalo de error. Suponer que los datos vienen afectados por un error máximo correspondiente al de aproximar por truncamiento las cifras dadas.

3. Aritmética de computadoresLos computadores no almacenan los números con precisión infinita sino de forma aproximada empleando un número fijo de bits (apócope del término inglés Binary Digit) o bytes (grupos de ocho bits). Prácticamente todos los computadores permiten al programador elegir entre varias representaciones o 'tipos de datos'. Los diferentes tipos de datos pueden diferir en el número de bits empleados, pero también (lo que es más importante) en cómo el número representado es almacenado: en formato fijo (también denominado 'entero') o en punto flotante2 (denominado 'real').

3.1 Aritmética de punto fijo

Un entero se puede representar empleando todos los bits de una palabra de computadora, con la salvedad de que se debe reservar un bit para el signo. Por ejemplo, en una máquina con longitud de palabra de 32 bits, los enteros están comprendidos entre -(231 - 1) y 231 - 1 = 2147483647. Un número representado en formato entero es 'exacto'. Las operaciones aritméticas entre números enteros son también 'exactas' siempre y cuando:

1.

La solución no esté fuera del rango del número entero más grande o más pequeño que se puede representar (generalmente con signo). En estos casos se dice que se comete un error de desbordamiento por exceso o por defecto (en inglés: Overflow y Underflow) y es necesario recurrir a técnicas de escalado para llevar a cabo las operaciones.

2.

La división se interpreta que da lugar a un número entero, despreciando cualquier resto.

Por estos motivos, la aritmética de punto fijo se emplea muy raramente en cálculos no triviales.

3.2 Números en punto flotante

.2.1 Notación científica normalizada

En el sistema decimal, cualquier número real puede expresarse mediante la denominada notación científica normalizada. Para expresar un número en notación científica normalizada multiplicamos o dividimos por 10 tantas veces como sea necesario para que todos los dígitos aparezcan a la derecha del punto decimal y de modo que el primer dígito después del punto no sea cero. Por ejemplo: 

En general, un número real x distinto de cero, se representa en notación científica normalizada en la forma: 

  (17)

en donde r es un número tal que  y n es un entero (positivo, negativo o cero).

Exactamente del mismo modo podemos utilizar la notación científica en el sistema binario. En este caso, tenemos que: 

  (18)

donde m es un entero. El número q se denomina mantisa y el entero m exponente. En un ordenador binario tanto q como m estarán representados como números en base 2. Puesto que la mantisa q está normalizada, en la representación binaria empleada se cumplirá que: 

 (19)

3.2.2 Representación de los números en punto flotante

En un ordenador típico los números en punto flotante se representan de la manera descrita en el apartado anterior, pero con ciertas restricciones sobre el número de dígitos de q y m impuestas por la longitud de palabra disponible (es decir, el número de bits que se van a emplear para almacenar un número). Para ilustrar este punto, consideraremos un ordenador hipotético que denominaremos MARC-32 y que dispone de una longitud de palabra de 32 bits (muy similar a la de muchos ordenadores actuales). Para representar un número en punto flotante en el MARC-32, los bits se acomodan del siguiente modo: 

Signo del número real x: 1 bit

Signo del exponente m: 1 bit

Exponente (entero |m|): 7 bits

Mantisa (número real |q|): 23 bits

 En la mayoría de los cálculos en punto flotante las mantisas se normalizan, es decir, se toman de forma que el bit más significativo (el primer bit) sea siempre '1'. Por lo tanto, la mantisa q cumple siempre la ecuación (19).

Dado que la mantisa siempre se representa normalizada, el primer bit en q es siempre 1, por lo que no es necesario almacenarlo proporcionando un bit significativo adicional. Esta forma de almacenar un número en punto flotante se conoce con el nombre de técnica del 'bit fantasma'.

Se dice que un número real expresado como aparece en la ecuación (18) y que satisface la ecuación (19) tiene la forma de punto flotante normalizado. Si además puede representarse exactamente con |m| ocupando 7 bits y |q| ocupando 24 bits, entonces es un número de máquina en el MARC-323

La restricción de que |m| no requiera más de 7 bits significa que:

 

Ya que  , la MARC-32 puede manejar números tan pequeños como 10-

38 y tan grandes como 1038. Este no es un intervalo de valores suficientemente generoso, por lo que en muchos casos debemos recurrir a programas escritos en aritmética de doble precisión e incluso de precisión extendida.

Como q debe representarse empleando no más de 24 bits significa que nuestros números de máquina tienen una precisión limitada cercana a las siete cifras decimales, ya que el bit menos significativo de la mantisa representa unidades de  . Por tanto, los números expresados mediante más de siete dígitos decimales serán objeto de aproximación cuando se almacenen en el ordenador.

Por ejemplo: 0.5 representado en punto flotante en el MARC-32 (longitud de palabra de 32 bits) se almacena en la memoria del siguiente modo:

 

Ejemplo 5: Suponga un ordenador cuya notación de punto fijo consiste en palabras de longitud 32 bits repartidas del siguiente modo: 1 bit para el signo, 15 bits para la parte entera y 16 bits para la parte fraccionaria. Represente los números 26.32,   y 12542.29301 en base 2 empleando esta notación de punto fijo y notación de punto flotante MARC-32 con 32 bits. Calcule el error de almacenamiento cometido en cada caso.

Solución: El número 26.32 en binario se escribe del siguiente modo:

Empleando las representaciones comentadas, obtenemos:

Si expresamos el error como la diferencia entre el valor y el número realmente almacenado en el ordenador, obtenemos:

En cuanto a los otros dos números, obtenemos:

 

Antes de entrar con detalle en la aritmética de los números en punto flotante, es interesante notar una propiedad de estos números de especial importancia en los cálculos numéricos y que hace referencia a su densidad en la línea real. Supongamos que p, el número de bits de la mantisa, sea 24. En el

intervalo  (exponente f = 0) es posible representar 224 números igualmente espaciados y separados por una distancia 1/224. De modo análogo, en cualquier

intervalo   hay 224 números equiespaciados, pero su densidad en este caso es 2f/224. Por ejemplo, entre 220 = 1048576 y 221 = 2097152 hay 224 =

16777216 números, pero el espaciado entre dos números sucesivos es de sólo  . De este hecho se deriva inmediatamente una regla práctica: cuando es necesario comparar dos números en punto flotante relativamente grandes, es siempre preferible comparar la diferencia relativa a la magnitud de los números. En la figura (1) se representa gráficamente la separación entre dos números consecutivos en función del exponente f en el rango f = [20,30].    

   

Figure: Evolución de la separación entre dos números consecutivos en función del exponente, f, de la representación en punto flotante de un número 

real.

[bb=55 60 455 410, clip=true, scale=0.7]eps/expon

3.3 Aritmética de punto flotante

En este apartado analizaremos los errores inherentes a la aritmética de los números de punto flotante. Primero consideraremos el error que surge como consecuencia de que los números reales no se pueden almacenar, en general, de forma exacta en un ordenador. Después analizaremos las cuatro operaciones aritméticas básicas y finalmente ampliaremos el estudio a un cálculo más complejo.

3.3.1 Números de máquina aproximados

Estamos interesados en estimar el error en que se incurre al aproximar un número real positivo x mediante un número de máquina del MARC-32. Si representamos el número mediante: 

en donde cada ai es 0 ó 1 y el bit principal es a1 = 1. Un número de máquina se puede obtener de dos formas:

Truncamiento: descartando todos los bits excedentes  . El número resultante, x' es siempre menor que x (se encuentra a la izquierda de x en la recta real).

Redondeo por exceso: Aumentamos en una unidad el último bit remanente a24 y después eliminamos el exceso de bits como en el caso anterior.

Todo lo anterior, aplicado al caso del MARC-32, se resume diciendo que si x es un número real distinto de 0 dentro del intervalo de la máquina, entonces el número de máquina x* más cercano a x satisface la desigualdad: 

 (20)

que se puede escribir de la siguiente forma: 

   

Ejemplo 6: ¿Cómo se expresa en binario el número x = 2/3? ¿Cuáles son los números de máquina x' y x'' próximos en el MARC-32?

El número 2/3 en binario se expresa como: 

Los dos números de máquina próximos, cada uno con 24 bits, son: 

x' =  

x'' =  

en donde x' se ha obtenido por truncamiento y x'' mediante redondeo por exceso. Calculamos ahora las diferencias x - x' y x'' - x para estimar cual es el error cometido: 

x - x' =  

x'' - x =  

Por tanto, el número más próximo es fl(x) = x'' y los errores de redondeo absoluto y relativo son: 

|fl(x) - x| =  

= 2-25 < 2-24  

3.3.2 Las operaciones básicas

Vamos a analizar el resultado de operar sobre dos números en punto flotante normalizado de l-dígitos de longitud, x e y, que producen un resultado normalizado del-dígitos. Expresaremos esta operación como: 

en donde op es +, -,   ó  . Supondremos que en cada caso la mantisa del resultado es primero normalizada y después redondeada (operación que puede dar lugar a un desbordamiento que requeriría renormalizar el número). El valor de la mantisa redondeada a p bits, qr, se define como (de una forma más rigurosa que en el caso anterior): 

en donde la función redondeo por defecto   es el mayor entero menor o igual 

a x y la función redondeo por exceso   es el menor entero mayor o igual a x. Para números enteros, esta función se traduce en la bien conocida regla de sumar 1 en la posición p + 1. Teniendo en cuenta sólo la mantisa, redondear de este modo da lugar a un intervalo máximo del error de: 

  (21)

y un error relativo máximo en el intervalo: 

 (22)

Analizaremos ahora el error generado por cada una de las operaciones básicas:

Multiplicación.

La operación de multiplicar dos números expresados en punto flotante implica sumar los exponentes y multiplicar las mantisas. Si la mantisa resultante no está normalizada, se recurre a renormalizar el resultado ajustando adecuadamente el exponente. Después, es necesario redondear la mantisa a p bits. Para analizar el error de esta operación supongamos dos números: 

Tenemos entonces que el producto será: 

xy = qx qy 2fx + f

y

en donde el valor de la mantisa se encontrará en el rango: 

ya que tanto x como y satisfacen la ecuación (19). Por tanto, la normalización del producto qx qy implica un desplazamiento a la derecha de, como máximo, una posición. La mantisa redondeada será entonces uno de estos dos posibles valores: 

en donde  , que es el error de redondeo, cumple la ecuación (21). Tenemos entonces: 

en donde, de acuerdo con la ecuación (22), tenemos: 

Por tanto, la cota del error relativo en la multiplicación es la misma que la que surge por redondear la mantisa.

División.

Para llevar a cabo la división en punto flotante, se divide la mitad de la mantisa del numerador por la mantisa del denominador (para evitar cocientes mayores de la unidad), mientras que los exponentes se restan. Esto es: 

Puesto que ambas mantisas satisfacen la ecuación (18), el valor del cociente estará acotado entre los límites: 

Aplicando un análisis similar al empleado en el caso de la multiplicación, obtenemos: 

en donde, de acuerdo con la ecuación (22), tenemos: 

Es decir, la cota máxima del error relativo en la división, como en el caso anterior, es la misma que la que surge por redondear la mantisa.

Adición y sustracción.

La operación de suma o resta se realiza del siguiente modo: se toma la mantisa del operando de menor magnitud (supongamos que es y) y se desplaza fx - fyposiciones a la derecha. La mantisa resultante es sumada (o restada) y el resultado se normaliza y después se redondea. Es decir: 

El análisis del error cometido en esta operación es más complejo que los estudiados hasta ahora, por lo que no lo vamos a ver en detalle. Sin embargo, el resultado final indica que la cota máxima del error cometido en la adición y la sustracción viene dado por: 

En conclusión, en todas las operaciones aritméticas elementales en punto flotante, el error absoluto del resultado es no mayor de 1 en el bit menos significativo de la mantisa.

Sin embargo, los errores de redondeo se acumulan a medida que aumenta el número de cálculos. Si en el proceso de calcular un valor se llevan a cabo Noperaciones aritméticas es posible obtener, en el mejor de los casos, un

error de redondeo total del orden de  4 (que coincide con el caso en que los errores de redondeo están aleatoriamente distribuidos, por lo que se produce una cancelación parcial). Desafortunadamente, este error puede crecer muy rápidamente por dos motivos:

Es muy frecuente que la regularidad del cálculo o las peculiaridades del computador den lugar a que el error se acumule preferentemente en una 

dirección; en cuyo caso el error de redondeo se puede aproximar a  . En circunstancias especialmente desfavorables pueden existir 

operaciones que incremente espectacularmente el error de redondeo. Generalmente, este fenómeno se da cuando se calcula la diferencia entre dos números muy próximos, dando lugar a un resultado en el cual los únicos bits significativos que no se cancelan son los de menor orden (en los únicos en que difieren). Puede parecer que la probabilidad de que se de dicha situación es pequeña, sin embargo, algunas expresiones matemáticas fomentan este fenómeno.

Veamos con un ejemplo los problemas comentados anteriormente. Hay dos formas de calcular las soluciones de la familiar ecuación cuadrática: 

ax2 + bx + c = 0

que son: 

 (23)

 (24)

Cualquiera de estas dos expresiones da problemas cuando a, c o ambos, son pequeños. En estos casos, el valor del discriminante es muy próximo al valor de b: 

por lo que la diferencia   viene afectada de un error de redondeo importante. En efecto, la ecuación (23) evalúa bien la raíz más grande en valor absoluto, pero da pésimos resultados al estimar la raíz menor en valor absoluto. Por otra parte, la ecuación (24) calcula bien la raíz menor (siempre en valor absoluto) pero no la raíz más grande.

La solución del problema pasa por emplear una expresión mejor condicionada. En este caso, es preferible calcular previamente: 

 (25)

y las dos raíces a partir de valor de q como: 

 (26)

Ejemplo: Calcular las raíces de la siguiente ecuación cuadrática: 

ax2 + bx + c = 0

siendo: 

Solución: Empleando la ecuación (23), obtenemos: 

Sin embargo, empleando la expresión (24): 

Por último, empleando las expresiones (25) y (26) se obtienen ambas soluciones correctas: 

3.4 Desbordamiento por exceso y desbordamiento por defecto

En la discusión anterior, hemos ignorado la posibilidad de que el resultado de una operación del punto flotante pueda no ser representable mediante el esquema fijo

(l-bits) empleado por el ordenador. La magnitud más grande que puede representarse mediante la fórmula general (18) es: 

  (27)

en donde F es el mayor exponente positivo representable (generalmente 27 - 1) y M es la mantisa que tiene todos sus bits puestos a 1 ( M = 1 - 224). Un desbordamiento por exceso de punto flotante (overflow en inglés) se origina cuando el resultado de una operación de punto flotante tiene una magnitud mayor que la representada por la ecuación (27).

Ejemplo: Con q = 8 (y por tanto F = 27 - 1 = 127), las siguientes operaciones aritméticas dan lugar a desbordamiento por exceso: 

El desbordamiento por defecto (underflow en inglés) se produce cuando el resultado de una operación en punto flotante es demasiado pequeño, aunque no nulo, como para que se pueda expresar en la forma dada por la ecuación (18). El número más pequeño representable suponiendo que siempre trabajamos con

mantisas normalizadas es  , en donde -F es el exponente negativo más grande permitido (generalmente -2-q-1). Por ejemplo, con q=8 resulta -F = -128.

Ejemplo: Con q = 8 (y por tanto -F = -128), la siguiente operación aritmética da lugar a desbordamiento por defecto: 

El desbordamiento por exceso es casi siempre resultado de un error en el cálculo. Sin embargo, en el caso del desobordamiento por defecto, en muchas ocasiones es posible continuar el cálculo reemplazando el resultado por cero.

3.5 Condicionamiento y estabilidadLa 'inestabilidad' en un cálculo es un fenómeno que se produce cuando los errores de redondeo individuales se propagan a través del cálculo incrementalmente. Veamos brevemente este fenómeno y el problema relacionado con este: el 'condicionamiento' del método o del problema.

La mejor forma de ver este fenómeno es a través de un ejemplo. Supongamos el siguiente sistema de ecuaciones diferenciales:

que tiene la siguiente solución general:

En el caso particular en que las condiciones iniciales de nuestro problema son:

y1(0) = -y2(0) = 1

es posible determinar que el valor de las constantes a1 y a2 es: a1 = 0 & y & a2 = 1

Hasta este punto, las soluciones son exactas. Sin embargo, supongamos que el sistema de ecuaciones anterior se resuelve empleando un método numérico cualquiera con el fin de calcular los valores de las funciones y1 y y2 en una

secuencia de puntos   y que el error del método da lugar a un valor

de  . Ya que a1 multiplica a un exponencial creciente cualquier valor, por pequeño que sea, de a1 dará lugar a que el término ex domine sobre el término e-

x para valores suficientemente grandes de x (ver figura (2)). La conclusión que se obtiene es que no es posible calcular una solución al sistema de ecuaciones diferenciales anterior que, para valores suficientemente grandes de x, no de lugar a un error arbitrariamente grande en relación con la solución exacta. 

   

   

Figure: Representación gráfica de las funciones y = e-

x e  en donde se pone de manifiesto que ambas funciones difieren rápidamente a partir de un cierto valor de la 

ordenada x.

[scale=0.7]eps/sinu

  

El problema anterior se dice que es inherentemente inestable, o empleando una terminología más común en cálculo numérico, se dice que está 'mal condicionado' (ill-conditioned).

4. Cálculo de raíces de ecuacionesEl objeto del cálculo de las raíces de una ecuación es determinar los valores de x para los que se cumple: 

 

f(x) = 0 (28)

La determinación de las raíces de una ecuación es uno de los problemas más antiguos en matemáticas y se han realizado un gran número de esfuerzos en este sentido. Su importancia radica en que si podemos determinar las raíces de una ecuación también podemos determinar máximos y mínimos, valores propios de matrices, resolver sistemas de ecuaciones lineales y diferenciales, etc...

La determinación de las soluciones de la ecuación (28) puede llegar a ser un problema muy difícil. Si f(x) es una función polinómica de grado 1 ó 2, conocemos expresiones simples que nos permitirán determinar sus raíces. Para polinomios de grado 3 ó 4 es necesario emplear métodos complejos y laboriosos. Sin embargo, si f(x) es de grado mayor de cuatro o bien no es polinómica, no hay ninguna fórmula conocida que permita determinar los ceros de la ecuación (excepto en casos muy particulares).

Existen una serie de reglas que pueden ayudar a determinar las raíces de una ecuación:

El teorema de Bolzano, que establece que si una función continua, f(x), toma en los extremos del intervalo [a,b] valores de signo opuesto, entonces la función admite, al menos, una raíz en dicho intervalo.

En el caso en que f(x) sea una función algebraica (polinómica) de grado n y coeficientes reales, podemos afirmar que tendrá n raíces reales o complejas.

La propiedad más importante que verifican las raíces racionales de una ecuación algebraica establece que si p/q es una raíz racional de la ecuación de coeficientes enteros: 

entonces el denominador q divide al coeficientes an y el numerador p divide al término independiente a0.

Ejemplo: Pretendemos calcular las raíces racionales de la ecuación: 

3x3 + 3x2 - x - 1 = 0

Primero es necesario efectuar un cambio de variable x = y/3: 

y después multiplicamos por 32: 

y3 + 3y2 -3y -9 = 0

con lo que los candidatos a raíz del polinomio son: 

Sustituyendo en la ecuación, obtenemos que la única raíz real es y = -3, es

decir,   (que es además la única raíz racional de la ecuación). Lógicamente, este método es muy poco potente, por lo que sólo nos puede servir a modo de orientación.

La mayoría de los métodos utilizados para el cálculo de las raíces de una ecuación son iterativos y se basan en modelos de aproximaciones sucesivas. Estos métodos trabajan del siguiente modo: a partir de una primera aproximación al valor de la raíz, determinamos una aproximación mejor aplicando una determinada regla de cálculo y así sucesivamente hasta que se determine el valor de la raíz con el grado de aproximación deseado.

.1 Método de la bisecciónEs el método más elemental y antiguo para determinar las raíces de una ecuación. Está basado directamente en el teorema de Bolzano explicado con anterioridad. Consiste en partir de un intervalo [x0,x1]tal que f(x0)f(x1) < 0, por lo que sabemos que existe, al menos, una raíz real. A partir de este punto se va reduciendo el intervalo sucesivamente hasta hacerlo tan pequeño como exija la 

precisión que hayamos decidido emplear.    

   

Figure: Diagrama de flujo correspondiente a la implementación del método de la bisección.

[scale=0.9]eps/bisecc

  

El algoritmo empleado se esquematiza en la figura (3). Inicialmente, es necesario suministrar al programa el número máximo de iteraciones MaxIter, la tolerancia  , que representa las cifras significativas con las que queremos obtener la solución y dos valores de la variable independiente, x0 y x1, tales que cumplan la relación f(x0)f(x1) < 0. Una vez que se comprueba que el intervalo de

partida es adecuado, lo dividimos en dos subintervalos tales que   

y   y determinamos en qué subintervalo se encuentra la raíz (comprobando de nuevo el producto de las funciones). Repetimos el proceso

hasta alcanzar la convergencia (hasta que  ) o bien hasta que se excede el número de iteraciones permitidas (Iter > MaxIter), en cuyo caso es necesario imprimir un mensaje de error indicando que el método no converge.

Dos operaciones representadas en el esquema de la figura (3) requieren una explicación adicional:

El punto medio del intervalo se calcula como   en lugar 

de emplear  . Se sigue de este modo una estrategia general al efectuar cálculos numéricos que indica que es mejor calcular una cantidad añadiendo un pequeño término de corrección a una aproximación obtenida previamente. Por ejemplo, en un computador de precisión limitada, existen valores de x0 y x1 para los cuales xm calculado 

mediante   se sale del intervalo [x0,x1]. La convergencia ( ) se calcula mediante la 

expresión  . De este modo, el término  , representa el número de cifras significativas con las que obtenemos el resultado.

4.2 Método de las aproximaciones sucesivasDada la ecuación f(x) = 0, el método de las aproximaciones sucesivas reemplaza esta ecuación por una equivalente, x=g(x), definida en la forma g(x)=f(x)+x. Para encontrar la solución, partimos de un valor inicial x0 y calculamos una nueva aproximación x1=g(x0). Reemplazamos el nuevo valor obtenido y repetimos el 

proceso. Esto da lugar a una sucesión de valores  , que si converge, tendrá como límite la solución del problema.    

  

Figure: Interpretación geométrica del método de las aproximaciones sucesivas.

[scale=0.9]eps/as-1

  

En la figura (4) se representa la interpretación geométrica del método. Partimos de un punto inicial x0 y calculamos y = g(x0). La intersección de esta solución con la recta y=x nos dará un nuevo valor x1 más próximo a la solución final.

Sin embargo, el método puede divergir fácilmente. Es fácil comprobar que el método sólo podrá converger si la derivada g'(x) es menor en valor absoluto que la unidad (que es la pendiente de la recta definida por y=x). Un ejemplo de este caso se muestra en la figura (5). Esta condición, que a priori puede considerarse una severa restricción del método, puede obviarse fácilmente. Para ello basta elegir la función g(x) del siguiente modo:

de forma que tomando un valor de   adecuado, siempre podemos hacer que g(x) cumpla la condición de la derivada.    

   

Figure: Demostración gráfica de que el método de las aproximaciones sucesivas diverge si la 

derivada g'(x) > 1.

[scale=0.9]eps/as-2

 4.3 Método de NewtonEste método parte de una aproximación inicial x0 y obtiene una aproximación mejor, x1, dada por la fórmula:

 (29)

 

La expresión anterior puede derivarse a partir de un desarrollo en serie de Taylor. Efectivamente, sea r un cero de f y sea x una aproximación a r tal que r=x+h. Sif'' existe y es continua, por el teorema de Taylor tenemos:

 

0 = f(r) = f(x+h) = f(x) + hf'(x) + O(h2)  (30)

 

en donde h=r-x. Si x está próximo a r (es decir hes pequeña), es razonable ignorar el término O(h2):

 

0 = f(x) + hf'(x)  (31)

 

por lo que obtenemos la siguiente expresión para h:

 (32)

 

A partir de la ecuación (32) y teniendo en cuenta que r=x+h es fácil derivar la ecuación (29).    

   

Figure: Interpretación geométrica del método de Newton.

[scale=0.9]eps/new-1

  

El método de Newton tiene una interpretación geométrica sencilla, como se puede apreciar del análisis de la figura (6). De hecho, el método de Newton consiste en una linealización de la función, es decir, f se reemplaza por una recta tal que contiene al punto (x0,f(x0)) y cuya pendiente coincide con la derivada de la función en el punto, f'(x0). La nueva aproximación a la raíz, x1, se obtiene de la intersección de la función linear con el eje X de ordenadas.

Veamos como podemos obtener la ecuación (29) a partir de lo dicho en el párrafo anterior. La ecuación de la recta que pasa por el punto (x0,f(x0)) y de pendientef'(x0) es:

 

y - f(x0) = f'(x0)(x-x0)  (33)

 

de donde, haciendo y=0 y despejando x obtenemos la ecuación de Newton-Raphson (29).    

   

Figure: Dos situaciones en las que el método de Newton no funciona adecuadamente: (a) el método no alcanza la convergencia y (b) el método converge hacia un punto que no es un cero de la ecuación.

[scale=0.9]eps/new-2

  

El método de Newton es muy rápido y eficiente ya que la convergencia es de tipo cuadrático (el número de cifras significativas se duplica en cada iteración). Sin embargo, la convergencia depende en gran medida de la forma que adopta la función en las proximidades del punto de iteración. En la figura (7) se muestran dos situaciones en las que este método no es capaz de alcanzar la convergencia (figura (7a)) o bien converge hacia un punto que no es un cero de la ecuación (figura (7b)).

4.4 Método de la secanteEl principal inconveniente del método de Newton estriba en que requiere conocer el valor de la primera derivada de la función en el punto. Sin embargo, la forma funcional de f(x) dificulta en ocasiones el cálculo de la derivada. En estos casos es más útil emplear el método de la secante.

El método de la secante parte de dos puntos (y no sólo uno como el método de Newton) y estima la tangente (es decir, la pendiente de la recta) por una aproximación de acuerdo con la expresión:

 (34)

 

Sustituyendo esta expresión en la ecuación (29) del método de Newton, obtenemos la expresión del método de la secante que nos proporciona el siguiente punto de iteración:

 (35)

      

   

Figure: Representación geométrica del método de la secante.

[scale=0.9]eps/secante

  

En la siguiente iteración, emplearemos los puntos x1 y x2para estimar un nuevo punto más próximo a la raíz de acuerdo con la ecuación (35). En la figura (8) se representa geométricamente este método.

En general, el método de la secante presenta las mismas ventajas y limitaciones que el método de Newton-Raphson explicado anteriormente.

.5 Método de Steffensen

El método de Steffensen presenta una convergencia rápida y no requiere, como en el caso del método de la secante, la evaluación de derivada alguna. Presenta además, la ventaja adicional de que el proceso de iteración sólo necesita un punto inicial. Este método calcula el siguiente punto de iteración a partir de la expresión: 

 (36)

4.6 Método de la falsa posiciónEl método de la falsa posición pretende conjugar la seguridad del método de la bisección con la rapidez del método de la secante. Este método, como en el método de la bisección, parte de dos puntos que rodean a la raíz f(x) = 0, es decir, dos puntos x0 y x1tales que f(x0)f(x1) < 0. La siguiente aproximación, x2, se calcula como la intersección con el eje X de la recta que une ambos puntos (empleando la ecuación (35) del método de la secante). La asignación del nuevo intervalo de búsqueda se realiza como en el método de la bisección: entre ambos intervalos, [x0,x2] y [x2,x1], se toma aquel que cumpla f(x)f(x2) < 0. En la figura (9) se representa geométricamente este método.    

   

Figure: Representación geométrica del método de la falsa posición.

[scale=0.9]eps/falpos

  

La elección guiada del intervalo representa una ventaja respecto al método de la secante ya que inhibe la posibilidad de una divergencia del método. Por otra parte y respecto al método de la bisección, mejora notablemente la elección del intervalo (ya que no se limita a partir el intervalo por la mitad).    

   

Figure: Modificación del método de la falsa posición propuesta por Hamming. La aproximación a la raíz se 

toma a partir del punto de intersección con el eje X de la recta que une los puntos ( x0,f(x0)/2) y 

(x1,f(x1)) si la función es convexa en el intervalo (figura a) o bien a partir de la recta que une los puntos 

(x0,f(x0)) y (x1, f(x1)/2) si la función es cóncava en el intervalo (figura b).

[scale=0.9]eps/hamming

  

Sin embargo, el método de la falsa posición tiene una convergencia muy lenta hacia la solución. Efectivamente, una vez iniciado el proceso iterativo, uno de los extremos del intervalo tiende a no modificarse (ver figura (9)). Para obviar este problema, se ha propuesto una modificación del método, denominada método de Hamming. Según este método, la aproximación a una raíz se encuentra a partir de la determinación del punto de intersección con el eje X de la recta que une los puntos ( x0,f(x0)/2) y (x1,f(x1)) si la función es convexa en el intervalo o bien a partir de la recta que une los puntos (x0,f(x0)) y (x1, f(x1)/2) si la función es cóncava en el intervalo. En la figura (10) se representa gráficamente el método de Hamming.

Como hemos comentado, el método de Hamming requiere determinar la concavidad o convexidad de la función en el intervalo de iteración. Un método relativamente sencillo para determinar la curvatura de la función consiste en evaluar la función en el punto medio del intervalo, f(xm) (en donde xm se calcula como en el método de la bisección) y comparar este valor con la media de los

valores de la función en los extremos del intervalo,  . Tenemos entonces que:

 5. Puntos fijos e iteración funcionalEl método de Newton y el de Steffenson son ejemplos de procedimientos mediante los cuales se calcula una sucesión de puntos empleando una fórmula de recurrencia como la siguiente: 

 

xn+1 = F(xn)   (37)

El algoritmo definido de este modo se denomina iteración funcional.

Ejemplo: En el caso del método de Newton, la expresión (37) se escribiría del modo: 

en tanto que para el método de Steffensen resulta ser: 

La fórmula (37) puede utilizarse para generar sucesiones que no convergen, como por ejemplo la sucesión 1, 3, 9, 27,..., que se obtiene con x0=1 y F(x)=3x. Sin embargo, estamos interesados en aquellos casos para los que

existe  . Es decir, aquellos casos para los que se cumple: 

Es fácil comprobar que si F es continua se cumple la siguiente relación entre F y s: 

Tenemos, por tanto, que F(s)=s y denominamos a s punto fijo de la función F. Podemos considerar al punto fijo como un valor al que se fija la función durante el proceso iterativo.

Como hemos visto en el apartado anterior, con frecuencia un problema matemático puede reducirse al problema de encontrar un punto fijo de una función. En este caso, nos limitaremos a analizar el caso más sencillo en

que F envía en sí mismo algún conjunto cerrado   y además se trata de una aplicación contractiva. Se dice que una transformación es contractiva si existe un número   menor que 1 que satisfaga la relación: 

  (38)

para todos los puntos x e y en el dominio de F.

Las aplicaciones contractivas cumplen una propiedad de gran importancia, que se puede expresar del siguiente modo:

Sea F una aplicación contractiva que va de un conjunto cerrado   a C. Entonces F tiene un punto fijo. Más aún, este punto fijo es el límite de toda sucesión que se obtenga a partir de la ecuación (37) con cualquier punto

inicial  .

El enunciado anterior, conocido como teorema de la aplicación contractiva se puede demostrar fácilmente. Para ello, primero escribimos xn en la forma: 

De acuerdo con la expresión anterior, vemos que la sucesión [xn]converge si y sólo si la serie 

converge. Para demostrar que esta serie converge, basta con demostrar que la serie 

 (39)

converge.

Por otra parte, usando la propiedad de las aplicaciones contractivas expresada por (38) junto con la ecuación (37), podemos escribir: 

  (40)

La relación expresada por (40) puede repetirse para obtener: 

  (41)

Para comprobar que la sucesión (39) converge, podemos utilizar el criterio de comparación, de modo que a partir de la expresión (41) obtenemos: 

Es decir, la sucesión converge tal como establece el teorema de la aplicación contractiva expresado anteriormente.

Comprobemos ahora que el punto fijo es efectivamente único. Para ello, supongamos que existen dos puntos fijos, x e y. De acuerdo con la relación (38), tenemos: 

Ya que   es un número finito menor que uno, la única forma de que la ecuación anterior se cumpla es si |x-y| = 0; es decir, si el punto fijo es único.

Ejemplo: Demuestre que la sucesión [xn] definida recursivamente de acuerdo con: 

es contractiva y tiene un punto fijo.

Para comprobar que la función anterior es contractiva, calculemos la diferencia entre dos términos cualesquiera de la sucesión anterior: 

(por la desigualdad triangular). Por tanto, de acuerdo con el teorema de la aplicación contractiva, la sucesión debe converger a un único punto fijo, cuyo valor es 2 (compruébelo).

6. Resolución de sistemas de ecuaciones linealesEl objetivo de este apartado es examinar los aspectos numéricos que se presentan al resolver sistemas de ecuaciones lineales de la forma: 

 

(42)

Se trata de un sistema de n ecuaciones con n incógnitas, x1, x2, ..., xn. Los elementos aij y bi son números reales fijados.

El sistema de ecuaciones (42) se puede escribir, empleando una muy útil representación matricial, como: 

 

(43)

Entonces podemos denotar estas matrices por A, x y b de forma que la ecuación se reduce simplemente a: 

 Ax=b (44)

Los métodos de resolución de sistemas de ecuaciones se pueden dividir en dos grandes grupos:

Los Métodos exactos o algoritmos finitos que permiten obtener la solución del sistema de manera directa.

Los Métodos aproximados que utilizan algoritmos iterativos e infinitos y que calculan las solución del sistema por aproximaciones sucesivas.

Al contrario de lo que pueda parecer, en muchas ocasiones los métodos aproximados permiten obtener un grado de exactitud superior al que se puede obtener empleando los denominados métodos exactos, debido fundamentalmente a los errores de truncamiento que se producen en el proceso.

De entre los métodos exactos analizaremos el método de Gauss y una modificación de éste denominado método de Gauss-Jordan. Entre los métodos aproximados nos centraremos en el estudio de los métodos de Richardson, Jacobi y Gauss-Seidel.

6.1 Métodos de resolución exacta

Antes de abordar el estudio de los métodos de resolución exacta de sistemas de ecuaciones lineales, analizaremos algunas propiedades y relaciones útiles que caracterizan a estos sistemas.

6.1.1 Sistemas fáciles de resolver

Analizaremos previamente un sistema que sea fácil de resolver. Por ejemplo,

supongamos que la matriz A de  presenta estructura diagonal, es decir, todos los componentes distintos de cero se encuentran sobre la diagonal principal. El sistema de ecuaciones (43) toma por tanto la forma: 

 

(45)

En este caso el sistema se reduce a n ecuaciones simples y la solución es: 

 

(46)

Continuando con la búsqueda de sistemas con soluciones fáciles, supongamos ahora que A tiene una estructura triangular inferior, es decir, todos los elementos de Adistintos de cero se sitúan bajo la diagonal principal: 

 

(47)

Es fácil ver que el valor de x1 se obtiene directamente a partir de la primera ecuación. Sustituyendo el valor conocido de x1 en la segunda ecuación es posible obtener el valor de x2. Procediendo de la misma forma para el resto de las ecuaciones, es posible obtener todos los valores x1 , x2, x3, ..., xn uno tras otro y en ese orden. El algoritmo formal para encontrar la solución se denomina sustitución progresiva y se puede expresar como: 

 

  (48)

Se puede emplear el mismo razonamiento para el caso en que la estructura de la matriz A sea triangular superior. En este caso el sistema matricial adopta la forma: 

 

(49)

y es posible obtener las soluciones en el orden xn, xn-1, ..., x1, empleando en este caso una modificación del algoritmo expresado por la ecuación (48) y que denominados algoritmo de sustitución regresiva: 

 

  (50)

Como es lógico, los métodos descritos se pueden aplicar a todos aquellos sistemas que se pueden convertir en un sistema triangular permutando filas y columnas de forma adecuada.

6.1.2 La factorización LU

Supongamos que A se puede factorizar como el producto de una matriz triangular inferior L con una matriz triangular superior U: 

 

A = LU (51)

En este caso, el sistema de ecuaciones dado por (44) podría representarse en la forma: 

 

LUx=b (52)

Si denominamos z a la matriz columna de n filas resultado del producto de las matrices Ux, tenemos que la ecuación (52) se puede reescribir del siguiente modo: 

 

Lz=b (53)

A partir de las ecuaciones (52) y (53), es posible plantear un algoritmo para resolver el sistema de ecuaciones empleando dos etapas:

Primero obtenemos z aplicando el algoritmo de sustitución progresiva en la ecuación (53).

Posteriormente obtenemos los valores de x aplicando el algoritmo de sustitución regresiva a la ecuación 

Ux = z

El análisis anterior nos muestra lo fácil que es resolver estos dos sistemas de ecuaciones triangulares y lo útil que resultaría disponer de un método que nos permitiera llevar a cabo la factorización A=LU. Si disponemos de una

matriz A de  , estamos interesados en encontrar aquellas matrices: 

tales que cumplan la ecuación (51). Cuando esto es posible, decimos que A tiene una descomposición LU. Se puede ver que las ecuación anterior no determina de forma única a Ly a U. De hecho, para cada i podemos asignar un valor distinto de cero a lii o uii (aunque no ambos). Por ejemplo, una elección simple es 

fijar lii=1 para   haciendo de esto modo que L sea una matriz triangular inferior unitaria. Otra elección es hacer U una matriz triangular superior unitaria (tomando uii=1 para cada i).

Para deducir un algoritmo que nos permita la factorización LU de Apartiremos de la fórmula para la multiplicación de matrices: 

 

(54)

en donde nos hemos valido del hecho de que lis=0 para s >i y usj=0 para s>j.

En este proceso, cada paso determina una nueva fila de U y una nueva columna

de L. En el paso k, podemos suponer que ya se calcularon las filas  

de U, al igual que las columnas   de L. Haciendo i=j=k en la ecuación (54) obtenemos 

 

(55)

Si especificamos un valor para lkk (o para ukk), a partir de la ecuación (55) es posible determinar un valor para el otro término. Conocidas ukk y lkk y a partir de la ecuación (54) podemos escribir las expresiones para la k-ésima fila (i=k) y para la k-ésima columna (j=k), respectivamente: 

 

  (56)

  (57)

Es decir, las ecuaciones (57) se pueden emplear para encontrar los elementos ukj y lik.

El algoritmo basado en el análisis anterior se denomina factorización de

Doolittle cuando se toman los términos lii = 1 para   (L triangular inferior unitaria) y factorización de Crout cuando se toman los términos uii=1 (U triangular superior unitaria).

Una implementación en pseudocódigo del algoritmo para llevar a cabo la factorización LU se muestra en la figura (11).

  

Figure: Implementación del algoritmo de la factorización LU.

Es interesante notar que los bucles que permiten el cómputo de la k-ésima fila de U y de la k-ésima columna de L se pueden llevar a cabo en paralelo, es decir, pueden evaluarse simultáneamente sobre dos procesadores, lo que redunda en un importante ahorro del tiempo de cálculo.

Ejemplo: Encuentre las factorizaciones de Doolittle y Crout de la matriz: 

La factorización de Doolittle es, a partir del algoritmo: 

En vez de calcular la factorización de Crout directamente, la podemos obtener a partir de la factorización de Doolittle que acabamos de ver. Efectivamente, si tenemos en cuenta que la matriz A es simétrica, es posible comprobar que se cumple la relación: 

A = LU = UTLT

por lo que la factorización de Crout resulta ser: 

6.1.3 Eliminación gaussiana básica

Ilustraremos el método de Gauss aplicando el procedimiento a un sistema de cuatro ecuaciones con cuatro incógnitas: 

 

(58)

En el primer paso, multiplicamos la primera ecuación por   y la restamos a 

la segunda, después multiplicamos la primera ecuación por   y la restamos 

a la tercera y finalmente multiplicamos la primera ecuación por   y la 

restamos a la cuarta. Los números 2,   y -1 son los multiplicadores del primer paso del proceso de eliminación. El número 6 es el elemento pivote de este 

primer paso y la primera fila, que no sufre modificación alguna, se denomina fila pivote. El sistema en estos momentos tiene el siguiente aspecto: 

 

(59)

En el siguiente paso del proceso, la segunda fila se emplea como fila pivote y -4 como elemento pivote. Aplicamos del nuevo el proceso: multiplicamos la 

segunda fila por   y la restamos de la tercera y después multiplicamos la 

segunda fila por   y la restamos a la cuarta. Los multiplicadores son en 

esta ocasión 3 y   y el sistema de ecuaciones se reduce a: 

 

(60)

El último paso consiste en multiplicar la tercera ecuación por   y restarla a la cuarta. El sistema resultante resulta ser: 

 

(61)

El sistema resultante es triangular superior y equivalente al sistema original (las soluciones de ambos sistemas coinciden). Sin embargo, este sistema es fácilmente resoluble aplicando el algoritmo de sustitución regresiva explicado en el apartado 6.1.1. La solución del sistema de ecuaciones resulta ser: 

Si colocamos los multiplicadores utilizados al transformar el sistema en una matriz triangular inferior unitaria (L) ocupando cada uno de ellos la posición del cero que contribuyó a producir, obtenemos la siguiente matriz: 

Por otra parte, la matriz triangular superior (U) formada por los coeficientes resultantes tras aplicar el algoritmo de Gauss (ecuación 61), es: 

Estas dos matrices nos dan la factorización LU de la matriz inicial de coeficientes, A, expresada por la ecuación (58): 

  

Figure: Implementación del algoritmo de eliminación gaussiana.

En la figura (12) se muestra un algoritmo en pseudocódigo para llevar a la práctica el proceso básico de eliminación gaussiana que acabamos de describir. En este algoritmo se supone que todos los elementos pivote son distintos de cero.

6.1.4 Método de Gauss-Jordan

Como hemos visto, el método de Gauss transforma la matriz de coeficientes en una matriz triangular superior. El método de Gauss-Jordan continúa el proceso de transformación hasta obtener una matriz diagonal unitaria (aij=0 para

cualquier  ).

Veamos el método de Gauss-Jordan siguiendo con el ejemplo empleado en el apartado anterior. Aplicando el método de Gauss habíamos llegado a la siguiente ecuación: 

Ahora seguiremos un procedimiento similar al empleado en el método de Gauss. Tomaremos como pivote el elemento a44=-3; multiplicamos la cuarta ecuación

por  y la restamos a la primera: 

Realizamos la misma operación con la segunda y tercera fila, obteniendo: 

Ahora tomamos como pivote el elemento a33=2, multiplicamos la tercera 

ecuación por   y la restamos a la primera: 

Repetimos la operación con la segunda fila: 

Finalmente, tomamos como pivote a22=-4, multiplicamos la segunda ecuación 

por   y la sumamos a la primera: 

El sistema de ecuaciones anterior es, como hemos visto, fácil de resolver. Empleando la ecuación (46) obtenemos las soluciones: 

6.1.5 Pivoteo

Sin embargo, los algoritmos de Gauss y Gauss-Jordan que acabamos de describir pueden dar lugar a resultados erróneos fácilmente. Por ejemplo, analicemos el

siguiente sistema de ecuaciones, en el que   es un número muy pequeño pero distinto de cero: 

Al aplicar el algoritmo gaussiano se obtiene el siguiente sistema triangular superior: 

y la solución es: 

En el computador, si   es suficientemente pequeño, los términos   

y   se computarán como un mismo número, por lo que   y  . Sin embargo, la solución correcta es: 

Tenemos entonces que la solución calculada es exacta para x2 pero extremadamente inexacta para x1.

El problema anterior no radica en la pequeñez del término aii, sino en su pequeñez relativa respecto de los otros elementos de su fila. La conclusión que podemos extraer es que un buen algoritmo debe incluir el intercambio de ecuaciones cuando las circunstancias así lo exijan. Un algoritmo que cumple este requisito es el denominado eliminación gaussiana con pivoteo de filas escaladas.

6.2 Métodos iterativos

El método de Gauss y sus variantes se conocen con el nombre de métodos directos: se ejecutan a través de un número finito de pasos y dan lugar a una solución que sería exacta si no fuese por los errores de redondeo.

Por contra, un método indirecto da lugar a una sucesión de vectores que idealmente converge a la solución. El cálculo se detiene cuando se cuenta con una solución aproximada con cierto grado de precisión especificado de antemano o después de cierto número de iteraciones. Los métodos indirectos son casi siempre iterativos: para obtener la sucesión mencionada se utiliza repetidamente un proceso sencillo.

6.2.1 Conceptos básicos

En general, en todos los procesos iterativos para resolver el sistema Ax=b se recurre a una cierta matriz Q, llamada matriz descomposición, escogida de tal forma que el problema original adopte la forma equivalente: 

 

Qx = (Q-A)x+b (62)

La ecuación (62) sugiere un proceso iterativo que se concreta al escribir: 

 (63)

El vector inicial x(0) puede ser arbitrario, aunque si se dispone de un buen candidato como solución éste es el que se debe emplear. La aproximación inicial que se adopta, a no ser que se disponga de una mejor, es la idénticamente 

nula  . A partir de la ecuación (63) se puede calcular una sucesión de vectores x(1), x(2), .... Nuestro objetivo es escoger una matriz Q de manera que:

se pueda calcular fácilmente la sucesión [x(k)]. la sucesión [x(k)] converja rápidamente a la solución.

Como en todo método iterativo, deberemos especificar un criterio de convergencia   y un número máximo de iteraciones M, para asegurar que el proceso se detiene si no se alcanza la convergencia. En este caso, puesto que x es un vector, emplearemos dos criterios de convergencia que se deberán satisfacer simultáneamente:

1.

El módulo del vector diferencia,  , partido por el módulo 

del vector x,   deberá ser menor que la convergencia deseada: 

2.

La diferencia relativa del mayor elemento en valor absoluto del 

vector x(k),  , deberá ser diez veces menor que  : 

6.2.2 Método de Richardson

El método de Richardson toma como matriz Q la matriz identidad (I). En este caso la ecuación (63) queda en la forma: 

 

Ix(k) = (I-A)x(k-1)+b = x(k-1)+r(k-1) (64)

en donde r(k-1) es el vector residual definido mediante r(k-1)=b-Ax(k-1).

La matriz identidad es aquella matriz diagonal cuyos elementos no nulos son 1, es decir: 

y cumple que 

IA = A

para cualquier valor de A; es decir, es el elemento neutro del producto matricial. De acuerdo con esto, la ecuación (64) se puede escribir como: 

x(k) = x(k-1) - Ax(k-1) + b = x(k-1) + r(k-1)

en donde un elemento cualquiera del vector r(k-1) vendrá dado por la expresión: 

En la figura (13) se muestra un algoritmo para ejecutar la iteración de Richardson. Este método recibe también el nombre de método de relajación o método de los residuos.

  

Figure: Implementación del algoritmo iterativo de Richardson.

6.2.3 Método de Jacobi

En la iteración de Jacobi, se escoge una matriz Q que es diagonal y cuyos elementos diagonales son los mismos que los de la matriz A. La matriz Q toma la forma: 

y la ecuación general (63) se puede escribir como 

 

Qx(k) = (Q-A)x(k-1) + b (65)

Si denominamos R a la matriz A-Q: 

la ecuación (65) se puede reescribir como: 

Qx(k) = -Rx(k-1) + b

El producto de la matriz Q por el vector columna x(k) será un vector columna. De modo análogo, el producto de la matriz R por el vector columna x(k-1) será también un vector columna. La expresión anterior, que es una ecuación vectorial, se puede expresar por necuaciones escalares (una para cada componente del vector). De este modo, podemos escribir, para un elemento i cualquiera y teniendo en cuenta que se trata de un producto matriz-vector: 

Si tenemos en cuenta que en la matriz Q todos los elementos fuera de la diagonal son cero, en el primer miembro el único término no nulo del sumatorio 

es el que contiene el elemento diagonal qii, que es precisamente aii. Más aún, los elementos de la diagonal de Rson cero, por lo que podemos eliminar el término i=j en el sumatorio del segundo miembro. De acuerdo con lo dicho, la expresión anterior se puede reescribir como:

de donde despejando xi(k) obtenemos: 

que es la expresión que nos proporciona las nuevas componentes del vector x(k) en función de vector anterior x(k-1) en la iteración de Jacobi. En la figura (14) se presenta un algoritmo para el método de Jacobi.

  

Figure: Implementación del método de Jacobi.

El método de Jacobi se basa en escribir el sistema de ecuaciones en la forma: 

 

(66)

Partimos de una aproximación inicial para las soluciones al sistema de ecuaciones y sustituimos estos valores en la ecuación (66). De esta forma, se genera una nueva aproximación a la solución del sistema, que en determinadas condiciones, es mejor que la aproximación inicial. Esta nueva aproximación se puede sustituir de nuevo en la parte derecha de la ecuación (66) y así sucesivamente hasta obtener la convergencia.

6.2.4 Método de Gauss-Seidel

La iteración de Gauss-Seidel se define al tomar Q como la parte triangular inferior de A incluyendo los elementos de la diagonal: 

Si, como en el caso anterior, definimos la matriz R=A-Q 

y la ecuación (63) se puede escribir en la forma: 

Qx(k) = -Rx(k-1) + b

Un elemento cualquiera, i, del vector Qx(k) vendrá dado por la ecuación: 

Si tenemos en cuenta la peculiar forma de las matrices Q y R, resulta que todos los sumandos para los que j > i en la parte izquierda son nulos, mientras que en 

la parte derecha son nulos todos los sumandos para los que  . Podemos escribir entonces: 

=  

=  

de donde despejando xi(k), obtenemos: 

Obsérvese que en el método de Gauss-Seidel los valores actualizados de xi sustituyen de inmediato a los valores anteriores, mientras que en el método de Jacobi todas las componentes nuevas del vector se calculan antes de llevar a cabo la sustitución. Por contra, en el método de Gauss-Seidel los cálculos deben llevarse a cabo por orden, ya que el nuevo valor xi depende de los valores actualizados de x1, x2, ..., xi-1.

En la figura (15) se incluye un algoritmo para la iteración de Gauss-Seidel.

  

Figure: Algoritmo para la iteración de Gauss-Seidel.

7. InterpolaciónNos centraremos ahora en el problema de obtener, a partir de una tabla de parejas (x,f(x)) definida en un cierto intervalo [a,b], el valor de la función para cualquierxperteneciente a dicho intervalo.

Supongamos que disponemos de las siguientes parejas de datos:

x x0 x1 x2 xn

y y0 y1 y2 yn

El objetivo es encontrar una función continua lo más sencilla posible tal que  

f(xi) = yi  (67)

Se dice entonces que la función f(x) definida por la ecuación (67) es una función de interpolación de los datos representados en la tabla.

Existen muchas formas de definir las funciones de interpolación, lo que da origen a un gran número de métodos (polinomios de interpolación de Newton, interpolación de Lagrange, interpolación de Hermite, etc). Sin embargo, nos centraremos exclusivamente en dos funciones de interpolación:

1.

Los polinomios de interpolación de Lagrange.2.

Las funciones de interpolación splines. Estas funciones son especialmente importantes debido a su idoneidad en los cálculos realizados con ordenador.

7.1 Polinomios de interpolación de Lagrange

Un polinomio de interpolación de Lagrange, p, se define en la forma: 

 

(68)

en donde   son polinomios que dependen sólo de los nodos 

tabulados  , pero no de las ordenadas  . La fórmula 

general del polinomio   es: 

 

(69)

Para el conjunto de nodos  , estos polinomios son conocidos como funciones cardinales. Utilizando estos polinomios en la ecuación (68) obtenemos la forma exacta del polinomio de interpolación de Lagrange.

Ejemplo: Suponga la siguiente tabla de datos:

x 5 -7 -6 0

y 1 -23 -54 -954

Construya las funciones cardinales para el conjunto de nodos dado y el polinomio de interpolación de Lagrange correspondiente.

Las funciones cardinales, empleando la expresión (69), resultan ser: 

El polinomio de interpolación de Lagrange es: 

7.2 Interpolación de splinesUna función spline está formada por varios polinomios, cada uno definido sobre un subintervalo, que se unen entre sí obedeciendo a ciertas condiciones de continuidad.

Supongamos que disponemos de n+1 puntos, a los que denominaremos nudos,

tales que  . Supongamos además que se ha fijado un

entero  . Decimos entonces que una función spline de grado k con nudos

en   es una función S que satisface las condiciones:

(i)

en cada intervalo  , S es un polinomio de grado menor o igual a k.

(ii)

S tiene una derivada de orden (k-1) continua en  .

Los splines de grado 0 son funciones constantes por zonas. Una forma explícita de presentar un spline de grado 0 es la siguiente:

Los intervalos   no se intersectan entre sí, por lo que no hay ambigüedad en la definición de la función en los nudos. Un spline de grado 1 se puede definir por:

 

En las figuras (16) y (17) se muestran las gráficas correspondientes a los splines de grado cero y de grado 1 respectivamente.    

   

Figure: Spline de grado 0 con seis puntos.

[scale=1.0]eps/spline-1

  

   

Figure: Spline de grado 1 con seis puntos.

[scale=1.0]eps/spline-2

 7.3 Splines cúbicosEl spline cúbico (k=3) es el spline más empleado, debido a que proporciona un excelente ajuste a los puntos tabulados y su cálculo no es excesivamente complejo.

Sobre cada intervalo  , S está definido por un polinomio cúbico diferente. Sea Si el polinomio cúbico que representa a S en el intervalo [ti,ti+1], por tanto:

Los polinomios Si-1 y Si interpolan el mismo valor en el punto ti, es decir, se cumple:

Si-1(ti) = yi = Si(ti)   

 

por lo que se garantiza que S es continuo en todo el intervalo. Además, se supone que S' y S'' son continuas, condición que se emplea en la deducción de una expresión para la función del spline cúbico.

Aplicando las condiciones de continuidad del spline S y de las derivadas primera S' y segunda S'', es posible encontrar la expresión analítica del spline. No

vamos a obtener esta expresión, ya que su demostración queda fuera del ámbito de estos apuntes. Simplemente diremos que la expresión resultante es:

En la expresión anterior, hi=ti+1-ti y   son incógnitas. Para determinar sus valores, utilizamos las condiciones de continuidad que deben cumplir estas funciones. El resultado (que tampoco vamos a demostrar) es:

La ecuación anterior, con   genera un sistema de n-1ecuaciones

lineales con n+1 incógnitas  . Podemos elegir z0 y z1 de forma arbitraria y resolver el sistema de ecuaciones resultante para obtener los valores

de  . Una elección especialmente adecuada es hacer z0=z1=0. La función spline resultante se denomina spline cúbico natural y el sistema de ecuaciones lineal expresado en forma matricial es:

 

(70)

 

en donde:

  

hi = ti+1-ti  

ui =  

bi =  

= (71)

      

  

Figure: Algoritmo para encontrar los coeficientes zi de un spline cúbico.

  

Este sistema de ecuaciones, que es tridiagonal, se puede resolver mediante eliminación gaussiana sin pivoteo como se muestra en la figura (18). El código acepta como entrada un conjunto de nodos (ti) y el conjunto de los valores de la función correspondiente (yi) y produce un vector con los vectores zi. Por último, el valor del spline S en un punto xcualquiera interpolado se puede calcular de forma eficiente empleando la siguiente expresión:

  (72)

 

en donde

Ai =  

Bi =  

Ci = (73)

  

Veamos un ejemplo para ilustrar el empleo de los splines cúbicos para interpolar los valores de una tabla. En la tabla (1) se muestran algunos valores de una serie

de valores tabulados a intervalos regulares de la función   en el intervalo [0,2.25]. También se indican los valores interpolados empleando el correspondiente spline cúbico así como el error absoluto cometido. Obsérvese que el error es cero para los nudos. En la figura (19) se representan gráficamente los valores tabulados.    

 

Table: Valores interpolados mediante un spline cúbico para la función   e indicación del error cometido (en valor absoluto).

x Si(x)

0.0000 0.0000 0.0000 0.0000E+00

0.0625   0.1426 1.0732E-01

0.1250   0.2782 7.5266E-02

0.1875   0.3997 3.3261E-02

0.2500 0.5000 0.5000 0.0000E+00

0.3125   0.5744 1.5440E-02

0.3750   0.6285 1.6155E-02

0.4375   0.6701 8.6732E-03

0.5000 0.7071 0.7071 0.0000E+00

       

1.7500 1.3228 1.3228 0.0000E+00

1.8125   1.3462 6.8994E-07

1.8750   1.3693 5.9953E-06

1.9375   1.3919 8.7004E-06

2.0000 1.4142 1.4142 0.0000E+00

2.0625   1.4361 2.4522E-05

2.1250   1.4577 4.7329E-05

2.1875   1.4790 4.6215E-05

2.2500 1.5000 1.5000 0.0000E+00

 

  

En la figura (20) se muestra otro ejemplo. Se representan gráficamente los puntos interpolados mediante una función spline cúbica para la función y=sen(x).    

   

Figure: Representación de la función  . Los círculos representan los valores tabulados de la 

función y la línea continua los puntos interpolados mediante una función spline cúbica.

[bb=560 195 70 580,scale=0.60,clip=true]eps/spline-3

  

  

Figure: Representación de la función y=sen(x). Los círculos representan los valores tabulados de la 

función y la línea continua los puntos interpolados mediante una función spline cúbica.

[bb=560 195 70 580,scale=0.60,clip=true]eps/spline-4

8. Integración numéricaDada una función f definida sobre un intervalo [a,b], estamos interesados en calcular 

 (74)

suponiendo que esta integral tenga sentido para la función f. La cuadratura o integración numérica consiste en obtener fórmulas aproximadas para calcular la integral J(f) de f. Estos métodos son de gran utilidad cuando la integral no se puede calcular por métodos analíticos, su cálculo resulta muy costoso y estamos interesados en una solución con precisión finita dada o bien sólo disponemos de una tabla de valores de la función (es decir, no conocemos la forma analítica de f).   8.1 Integración vía interpolación polinomial

Una estrategia muy útil para calcular el valor numérico de la integral dada por la ecuación (74) consiste en reemplazar fpor otra función g, fácil de integrar, que

aproxima a f de forma adecuada. Si  , se deduce que 

Los polinomios son buenos candidatos para el papel de g. De hecho, g puede ser un polinomio que interpola a f en cierto conjunto de nodos5.

Supongamos que deseamos calcular la integral (74). Podemos elegir una serie de

nudos,   en el intervalo [a,b] e iniciar un proceso de interpolación de Lagrange (ver apartado 7.1 para una descripción de los polinomios de interpolación de Lagrange). El polinomio de grado menor o igual a n que interpola a f en los nudos es: 

(75)

La integral (74) se puede escribir entonces como: 

Es decir, tenemos una fórmula general que se puede emplear para cualquier f y que tiene la forma: 

 (76)

en donde 

8.2 Regla del trapecio

Si en la expresión (76) empleamos polinomios de grado n=1 y tomamos como nudos x0=a y x1=b, tenemos el caso más sencillo posible, en donde los polinomios de interpolación son: 

=  

=  

por lo que: 

La fórmula de cuadratura correspondiente es: 

Esta expresión se conoce como regla del trapecio y proporciona un resultado exacto para todas las funciones de grado menor o igual a 1.

8.3 Regla de Simpson

Empleando un razonamiento similar al anterior y tomando un polinomio de grado n=2 para interpolar a f, obtenemos la conocida regla de Simpson: 

 (77)

que es exacta para todos los polinomios de grado   2 y curiosamente, exacta 

para todos los polinomios de grado   3.

En los cálculos prácticos se emplea, generalmente, la regla de Simpson compuesta, en la que el intervalo de integración [a,b] se divide en un número par, n, de subintervalos. Tenemos entonces: 

en donde 

h = (b-a)/n

Aplicando la regla de Simpson (77) en cada uno de los subintervalos se obtiene la expresión final: 

(78)