TRABAJO FINAL DE METODOS NUMERICOS

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 diversos1

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)

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

= (6)
 = (7)

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 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 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 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. 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.

3.1 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

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.

3.2 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:

(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.

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)).

3.3 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:
 (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.

3.4 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)

3.5 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.

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:

4. 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:

(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.

4.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.

4.1.1 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 lak-é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:

4.1.2 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:
4.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.

4.2.1 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.

4.2.2 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.3 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.

5. 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.

5.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 deLagrange.

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:

5.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:

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

6. Integración numérica
Dada 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).

6.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

6.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.

6.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)