Libro Métodos
numéricos Wladimiro
Diaz Villanueva
Se puede descargar de los links:
Transcrito
a Word por: Alexander Arias Londoño
Fecha: 3
de junio de 2017
- 1. Introducción
- 2. Errores
- 3. Aritmética de computadores
- 3.1 Aritmética de punto fijo
- 3.2 Números en punto flotante
- 3.3 Aritmética de punto flotante
- 3.4 Desbordamiento por exceso y desbordamiento por defecto
- 3.5 Condicionamiento y estabilidad
- 4. Cálculo de raíces de ecuaciones
- 4.1 Método de la bisección
- 4.2 Método de las aproximaciones sucesivas
- 4.3 Método de Newton
- 4.4 Método de la secante
- 4.5 Método de Steffensen
- 4.6 Método de la falsa posición
- 5. Puntos fijos e iteración funcional
- 6. Resolución de sistemas de ecuaciones lineales
- 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
- 6.2 Métodos iterativos
- 7. Interpolación
- 8. Integración numérica
- About this document ...
1. Introducción
La 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 diversos.
Un ejemplo: durante la segunda
guerra mundial, el Ballistic Research Laboratory (BRL), un laboratorio
del ejército americano, destinaba todos sus esfuerzos al cálculo de ``tablas de
tiro'', necesarias para ajustar el ángulo de tiro en función de numerosos factores
(tipo de proyectil y cañón, situación del blanco, velocidad del viento,
temperatura, etc...). Estas tablas eran largas y difíciles de establecer: sólo
con dos factores (alcance del proyectil y altitud del blanco) se requerían
entre 2000 y 4000 trayectorias, cada una de las cuales requería 750
multiplicaciones de 10 cifras. Una ``calculadora humana'' requería tres días
para calcular una sola trayectoria. En el verano de 1944 el BRL producía quince
tablas por semana mientras que la demanda era de cuarenta, por lo que se
paralizó la entrada en servicio de nuevos proyectiles, ya que de nada servían
en el frente sin las famosas tablas. En febrero de 1946 (tan sólo año y medio
después), el primer ordenador ``moderno'', denominado ENIAC, era capaz de calcular
una trayectoria cada veinte segundos.
2. Errores
El 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)
|
|
|
=
|
|
(6)
|
|
|
=
|
|
2.2 Dígitos significativos
Sea
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 d decimales, 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
|
|
|
|
=
|
|
|
|
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)
|
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:
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:
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 adicionales
1.
¿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 F
perpendicular 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 computadores
Los 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 flotante (denominado 'real').
La
aritmética de punto flotante que se emplea actualmente en todos los
ordenadores, fue puesta a punto por el ingeniero español Leonardo Torres y
Quevedo (1852-1936). Entre sus logros se encuentran, además de varios autómatas
teledirigidos y una máquina que jugaba ajedrez (y que nunca perdía), una
calculadora que se manejaba mediante una máquina de escribir. En 1914 presentó
los planos completos de una <<máquina analítica>>.
Desafortunadamente, no se llegó a construir jamás.
- 3.1 Aritmética de punto fijo
- 3.2 Números en punto flotante
- 3.3 Aritmética de punto flotante
- 3.4 Desbordamiento por exceso y desbordamiento por defecto
- 3.5 Condicionamiento y estabilidad
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
3.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:
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:
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:
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-32
La mayoría de los números reales no pueden
representarse de manera exacta en la MARC-32.
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:
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:
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 de l-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:
y un error
relativo máximo en el intervalo:
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 + fy
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:
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:
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
- fy posiciones 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:
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 N operaciones aritméticas es posible obtener, en el mejor de los casos, un error de redondeo total del orden de (que coincide con el caso en que los errores de redondeo están aleatoriamente distribuidos, por lo que se produce una cancelación parcial).
En realidad, es el número en punto flotante más pequeño que, cuando se suma a 1.0, da lugar a un número diferente de 1.0. Se denomina 'precisión de máquina' y a todos los efectos consideraremos que su valor es muy próximo al error de redondeo comentado. A efectos comparativos, si empleamos palabras de 32 bits su valor es de aproximadamente
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.
ax2 + bx + c = 0
que son:
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:
y las dos
raíces a partir de valor de q como:
ax2 + bx + c = 0
siendo:
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:
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: Ejemplo: Con q = 8 (y por tanto -F = -128), la siguiente operación aritmética da lugar a desbordamiento por defecto:
3.5 Condicionamiento y estabilidad
La
'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
|
4. Cálculo de raíces de ecuaciones
El 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.
- 4.1 Método de la bisección
- 4.2 Método de las aproximaciones sucesivas
- 4.3 Método de Newton
- 4.4 Método de la secante
- 4.5 Método de Steffensen
- 4.6 Método de la falsa posición
4.1 Método de la bisección
Es 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
|
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 sucesivas
Dada 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
|
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:
|
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 Newton
Este
método parte de una aproximación inicial x0 y obtiene una
aproximación mejor, x1, dada por la fórmula:
|
0 = f(r) = f(x+h)
= f(x) + hf'(x) + O(h2)
|
(30)
|
|
0 = f(x) + hf'(x)
|
(31)
|
|
Figure: Interpretación
geométrica del método de Newton.
|
|
[scale=0.9]eps/new-1
|
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 pendiente f'(x0) es:
|
y - f(x0)
= f'(x0)(x-x0)
|
(33)
|
|
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
|
4.4 Método de la secante
El
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: |
Figure: Representación
geométrica del método de la secante.
|
|
[scale=0.9]eps/secante
|
En general, el método de la secante presenta las mismas ventajas y limitaciones que el método de Newton-Raphson explicado anteriormente.
4.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:4.6 Método de la falsa posición
El 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
|
|
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
|
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 funcional
El 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:
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:
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
converge.
Por otra parte, usando la propiedad de las aplicaciones contractivas
expresada por (38)
junto con la ecuación (37), podemos
escribir:
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 lineales
El
objetivo de este apartado es examinar los aspectos numéricos que se presentan
al resolver sistemas de ecuaciones lineales de la forma:
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:
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
- 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
- 6.2 Métodos iterativos
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
- 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
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:
En este
caso el sistema se reduce a n ecuaciones simples y la solución es:
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)
|
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)
|
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)
|
|
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
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:
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
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)
|
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.
|
|
|
Ejemplo: Encuentre las factorizaciones de Doolittle y Crout de la matriz:
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:
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:
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:
El último
paso consiste en multiplicar la tercera ecuación por
y restarla a
la cuarta. El sistema resultante resulta ser:
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:
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.
|
|
|
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:
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
- 6.2.2 Método de Richardson
- 6.2.3 Método de Jacobi
- 6.2.4 Método de Gauss-Seidel
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)
|
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.
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:|
Qx(k) = (Q-A)x(k-1)
+ b
|
(65)
|
Si
denominamos R a la matriz A-Q:
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.
|
|
|
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
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:
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ón
Nos
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 cualquier xperteneciente 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:
en donde
son
polinomios que dependen sólo de los nodos tabulados
, pero no de
las ordenadas
. La fórmula
general del polinomio
es:
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: 7.2 Interpolación de splines
Una 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:
|
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úbicos
El 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: |
Si-1(ti) = yi
= Si(ti)
|
|
|
|
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:
|
hi
|
=
|
ti+1-ti
|
|
|
ui
|
=
|
|
|
|
bi
|
=
|
|
|
|
|
=
|
|
(71)
|
|
Figure: Algoritmo para
encontrar los coeficientes zi de un spline cúbico.
|
|
|
|
Ai
|
=
|
|
|
|
Bi
|
=
|
|
|
|
Ci
|
=
|
|
(73)
|
|
Table: Valores
interpolados mediante un spline cúbico para la función
e indicación del error cometido (en valor
absoluto).
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
No hay comentarios:
Publicar un comentario