La resolución de ecuaciones y sistemas (ya sean polinómicas, algebraicas o trascendentes) es uno de los problemas que con más frecuencia aparece en los distintos campos de la Ciencia y la Técnica. Además se trata de un problema que ha sido estudiado desde muy antiguamente; por ejemplo ya en el año 100 A.C. Herón empleaba un método iterativo para aproximar la raíz cuadrada de un número positivo.
No existen métodos generales de resolución simbólica de ecuaciones o sistemas; salvo para ciertas ecuaciones de tipo polinómico o en el caso de sistemas lineales. Por ello surge la necesidad de desa--rrollar métodos numéricos para calcular, al menos de forma aproximada, las soluciones de este tipo de problemas. Teniendo en cuenta estas cuestiones, Mathematica incorpora diferentes comandos para la re-solución de ecuaciones y sistemas de ecuaciones de modo formal (por ejemplo, Solve, Reduce, Eliminate, Roots), siempre que admitan solución simbólica; pero en la mayoría de los casos se calculará una aproximación numérica de dicha solución me-dian-te otras órdenes (por ejemplo, NSolve, NRoots, FindRoot), empleándose técnicas apropiadas para cada pro-blema concreto.
Comenzaremos describiendo alguno de los métodos numéricos más conocidos de
aproximación de soluciones para una ecuación ge-ne-ral, de la forma para
al menos continua y con
distinto signo en los extremos de
.
Continuaremos después con ecuaciones de tipo polinómico, tanto con la descripción de comandos directos de Mathematica para su resolución formal y aproximada, como con ciertos algoritmos y particularidades específicas para este tipo de ecuaciones.
Pretendemos en esta sección describir los algoritmos de algunos métodos numé-ricos de resolución aproximada de ecuaciones y estudiar las condiciones en las que éstos determinan una sucesión de aproximaciones convergente a la solución real.
Ya que, en general, no existen métodos directos (ni fórmulas, ni algoritmos)
para calcular formalmente las raíces de una ecuación del tipo (ecgen) y sólo un número muy reducido de ecuaciones
pueden resolverse de forma exacta (por ejemplo, la ecuación de segundo grado),
los métodos numéricos que se emplean para resolver estas ecuaciones son de tipo
iterativo y proporcionan una sucesión de aproximaciones, . Lo interesante será que esta
sucesión de aproximaciones converja hacia una raíz, pongamos
, de la
ecuación.
Vamos a estudiar alguno de los métodos más empleados en la práctica, como son el de bisección, el de la secante, regula-falsi, etc., y el de Newton-Raphson. El primero requiere menos hipótesis para asegurar la convergencia, pero en cambio es más lento que el último; es decir, si queremos obtener una aproximación con una precisión determinada, en general tendremos que hacer más iteraciones con el primero. Sin embargo, para poder aplicar el de Newton-Raphson no basta con que la función sea continua, como veremos después.
Es el método más simple que se puede emplear para resolver ecuaciones. Sólo
requiere que la función sea continua y que hayamos localizado un cambio de signo
de la misma en los extremos de cierto intervalo en el que empezaremos a
trabajar. Consiste, pues, en la aplicación reiterada del conocido teorema de
Bolzano, una vez asegurada la existencia de al menos una solución de la ecuación
(ecgen) en el intervalo .
El único inconveniente que surge es la posible existencia de más de una de estas
raíces, ya que el teorema de Bolzano no nos asegura la unicidad de la misma.
Pero en todo caso el correspondiente algoritmo siempre conducirá a la
aproximación numérica de una de estas raíces. Nos podemos ayudar de una
representación gráfica previa de la función para seleccionar un intervalo en el
que sólo exista una de estas raíces.
Recordamos el teorema fundamental en el que se basa dicho algoritmo
Sea una función continua en
y tal que
. Entonces existe
tal que
.
En la figura c02g1 se presenta la interpretación
geométrica. Se considera el intervalo . En los extremos la función
cambia de signo, y hay más de un cero.
Idea geométrica del teorema de
Bolzano.
Supongamos que verifica las hipótesis del
teorema de Bolzano en cierto subintervalo
del intervalo de
partida
.
Entonces, también podemos asegurar la existencia de al menos una raíz de la
ecuación,
.
Sea ahora , que es el punto medio
del intervalo. En valor absoluto, el error cometido al tomar
como
aproximación de
(ver figura c02g2 con
y
), es inferior a
la mitad de la longitud de dicho intervalo, es decir
. Lo indicamos
como sigue:
El cálculo del punto medio se
efectúa mediante la fórmula
pues es más conveniente
en términos de implementación en un ordenador, cuya precisión será siempre
limitada (ver la justificación y estrategia general en Kincaid-Cheney).
Como es
conocida la podemos evaluar en
y ver si este punto es una
raíz exacta de la ecuación, en cuyo caso ya habríamos finalizado la búsqueda;
o bien, en caso contrario, podemos determinar el signo de la función en el
centro de este intervalo, de tal forma que dicho signo será opuesto del que
tiene
en uno de
los dos extremos, ya que
y
Aquí hay que
hacer notar que, a no ser que trabajemos con un software que permita el
cálculo simbólico, debido a los inevitables errores de redondeo en la
representación interna de los números en coma flotante en la memoria del
ordenador y los acumulados mediante las sucesivas operaciones efectuadas, no
conviene implementar una condición del tipo ?`
? sino más bien
?`
para cierto valor
de
,
suficientemente pequeño?
Así, eligiendo el extremo en el que tiene signo opuesto que en
el centro, tenemos un nuevo intervalo que denotaremos
, de tamaño mitad
que el anterior, en el cual seguimos teniendo asegurado que
tiene una raíz. La nueva
aproximación será el centro de este nuevo intervalo y el error máximo cometido
es la mitad de la longitud del mismo, es decir, la cuarta parte de la longitud
del intervalo anterior. En otros términos,
En general, empezando para con el intervalo de
partida, y mediante un proceso de inducción, tras
iteraciones, el valor
absoluto del error cometido satisface la desigualdad
Vemos pues cómo, si no se llega a encontrar la solución exacta para algún
, entonces se
obtiene una sucesión de intervalos encajados
con
, de manera que, si
para todo
, entonces el error
tiende a cero cuando
tiende a infinito, y, por la
complitud del cuerpo de los números reales, el método
resulta convergente, es decir
en
. Como es
lógico, en la práctica computacional no se lleva a cabo dicho paso al límite,
pues implicaría realizar un número infinito de iteraciones (bucle
infinito); además, superado un determinado número de éstas, ya estaríamos
llegando incluso a entrar en conflicto con el grado de precisión de la máquina
con que estemos trabajando (consultar capítulo anterior sobre los errores de
redondeo). Mathematica permite controlar esto mediante
diferentes órdenes, como son WorkingPrecision, SetPrecision,
AccuracyGoal, etc. Así pues, lo que se suele hacer es intentar asegurar
una determinada precisión o tolerancia, que denominaremos
, en la aproximación; para
ello debemos detener el proceso cuando la longitud del intervalo correspondiente
sea menor o igual que
; también podemos calcular de
antemano el número de iteraciones resolviendo la inecuación
, lo que es inmediato tomando
logaritmos.
Describimos a continuación el algoritmo en forma de pseudocódigo. Aquí, con el afán de simplificar y optimizar al
máximo el código correspondiente y la cantidad de variables a almacenar lo que
se hace es mantener en todo momento las letras de las variables y
para los extremos del
intervalo considerado y
para el punto medio del
mismo. En la figura c02g2 se ilustra.
Idea geométrica del método de
bisección
ALGORITMO DE BISECCIÓN |
Entrar ![]() ![]() ![]() ![]() |
Mientras ![]() |
Hacer ![]() |
Si ![]() ![]() |
Si ![]() ![]() |
Si ![]() ![]() |
En el método de regula falsi se puede partir de las
mismas hipótesis que en el de bisección, suponiendo que se verifica el teorema
de Bolzano en el intervalo ; sólo que ahora el
algoritmo para ir aproximando a una de las raíces es ir tomando la
correspondiente recta secante a la gráfica en los puntos
y
,
, y tomar la
intersección de ésta con el eje
(ver figura c02g3); llamemos también
a la aproximación de la raíz
buscada. A continuación efectuaremos un chequeo de cambio de signo, lo mismo que
en el método de bisección, y nos quedaremos con el subintervalo donde se siga
manteniendo el cambio de signo,
o bien
. Bastará con
repetir este proceso de forma recursiva para ir obteniendo cada vez mejores
aproximaciones de la raíz buscada.
Idea geométrica del método de
regula-falsi.
Así pues, el correspondiente algoritmo para el método de regula-falsi va a
resultar muy parecido al del método de bisección, ya presentado anteriormente,
cambiando solamente el paso en el que se calcula como el punto de corte de la
correspondiente recta secante, de ecuación
con el eje
:
Se trata de una operación de
cálculo elemental que obviamente podemos realizar a mano sin ningún problema;
pero también es posible utilizar las facilidades de cálculo simbólico que poseen
los paquetes de software matemático actuales como Mathematica para realizar este tipo de desarrollos rutinarios
elementales.
En cuanto a la convergencia de este nuevo método, se siguen obteniendo
intervalos encajados (sólo que ahora su longitud no tiene que tender hacia cero)
y una sucesión de aproximaciones de manera que
, siendo
una raíz, que sabemos que
posee la ecuación (ecgen) en el intervalo de
partida.
Geométricamente este método se basa en la misma idea que el de regula-falsi,
salvo que ahora se hará caso omiso a la sucesión de intervalos encajados que
contienen a la raíz de la ecuación (ecgen) y
simplemente se seguirá un proceso iterativo a partir de los valores iniciales
y
mediante la
siguiente fórmula (obtenida de la misma forma que en (metrfalsi), con las identificaciones
y
):
De esta manera nos evitamos el tener que chequear en cada paso del algoritmo, el correspondiente cambio de signo de la función en los extremos de los intervalos, pero corremos el riesgo de que en algún caso no se tenga la deseada convergencia hacia la raíz de la ecuación (ecgen).
En el método de Whittaker, partiendo de un valor
inicial , se
aproxima la raíz de la ecuación (ecgen)
mediante la abscisa del correspondiente punto de corte con el eje
de la recta que pasa por los
sucesivos puntos
,
, y tiene como
pendiente un valor prefijado
fijo (ver figura c02g4 ).
Idea geométrica del método de
Whittaker
Así pues, la ecuación de esas rectas es y el cálculo de
a partir de
por
necesitándose en este caso hipótesis bastante más restrictivas que en los
métodos anteriores para poder asegurar la convergencia. Se pretende tomar un
valor
que
aproxime, para todas las iteraciones, el valor
de (metsec) y, por tanto, simplificar los cálculos de
cada iteración.
Este método se basa en la sustitución de la función en la ecuación (ecgen) por un polinomio de segundo grado que pase
por tres puntos dados,
,
y
(con
), de la
correspondiente gráfica
hallando posteriormente el
punto de corte de dicha parábola con el eje
mediante la conocida
fórmula para este tipo ecuaciones de segundo grado (ver figura c02g5 ).
Idea geométrica del método de
Müller.
Este método es menos usado pues, entre otras cosas, presenta el problema de elegir cuál de las dos raíces de la ecuación de segundo grado es la adecuada en cada paso.
En este caso se obtiene una sucesión de aproximaciones partiendo de un valor
inicial , que
debe ser dado o elegido convenientemente. Una vez conocida una aproximación,
digamos
, la
siguien-te,
, se
obtiene hallando el punto de corte de la correspondiente recta tangente a la
curva de ecuación
en el punto
con el eje
(ver figura c02g6 ).
Idea geométrica del método de
Newton-Raphson.
El proceso se repite sucesivamente, obteniéndose pues el siguiente método
iterativo: para dado (bastará con tomar
en (metWhittaker), con la salvedad de que ahora esta
pendiente irá cambiando en cada iteración)
El método iterativo empleado ya por Herón 100 años antes del comienzo de
nuestra era para aproximar la raíz cuadrada de un número positivo ,
es el método que acabamos
de deducir aplicado a
.
Además podemos observar que el método de la secante (metsec) puede obtenerse como una modificación del
método de Newton-Raphson, en el que el valor de la pendiente de la recta
tangente se aproxima por el valor de la pendiente de la correspondiente recta
secante
Para una curva que no tenga grandes oscilaciones alrededor de cierto punto,
la recta tangente puede resultar una buena aproximación de la misma, al menos en
un cierto entorno de dicho punto, y, por tanto, se puede ver que el método va a
converger en la mayoría de las ocasiones. En general, esto ocurre siempre que se
parta de una buena aproximación inicial (como se comprobará más
adelante). Una práctica habitual consiste, una vez detectado un cambio de signo
en la función, en emplear el método de bisección para obtener una primera
aproximación de la raíz de la ecuación, que es utilizada como valor inicial
del método de
Newton-Raphson .
Hemos de hacer algunas observaciones al método. En primer lugar, debe ser derivable para
poder trazar la recta tangente. Por otra parte, si la derivada se anula en algún
punto y dicho punto aparece como una de las sucesivas aproximaciones, no podemos
continuar, ya que en la siguiente aparecería una división por cero. E incluso
aunque no se divida entre cero, la división entre números muy próximos a cero
puede dar lugar a grandes errores de redondeo que afectarán en gran medida a la
precisión de los resultados.
Por último, el control del error no es tan inmediato como en el método de bisección. En este caso suele utilizarse como método de parada la distancia entre dos aproximaciones consecutivas.
Planteamos ya el algoritmo de Newton-Raphson. Con el mismo afán de
simplificación que antes, reutilizamos las variables y
para almacenar
los valores de
y
,
respectivamente.
ALGORITMO DE NEWTON-RAPHSON |
Entrar ![]() ![]() ![]() |
Hacer ![]() |
Si ![]() ![]() |
En caso contrario intercambiar ![]() ![]() |
A continuación damos un resultado sobre convergencia del método de Newton--Raphson, que es fácil de comprobar en muchos ejemplos.
Supongamos que la función admite derivada segunda en
y verifica las siguientes propiedades:
para todo
.
no cambia de signo en
(por ejemplo,
en todo
punto).
Entonces, partiendo de un valor inicial para
el que se
verifique que , la sucesión de
aproximaciones obtenida por el método de Newton-Raphson converge hacia la
única raíz de la ecuación
en
.
Antes de probar el teorema, observemos que bastará con tomar por ejemplo como
el extremo del
intervalo para el cual
y
tengan el mismo signo. Nótese, además, que la existencia de la raíz está
asegurada, gracias al teorema de Bolzano, por la continuidad de
y la primera hipótesis.
La unicidad de la raíz se debe a la segunda, ya que la función será estrictamente monótona (creciente o decreciente) en todo el intervalo. La tercera hipótesis garantiza que la recta tangente no corta a la curva, ya que ésta tiene concavidad fija.
En cuanto a las cuatro distintas posibilidades geométricas que se pueden dar
con estas hipótesis (ver figura c02g7), es inmediato
comprobar que todas pueden reducirse a una cualquiera de ellas sin más que tomar
la reflexión adecuada respecto al eje (mediante el cambio
) o al eje
(mediante el cambio
, que
transformaría a su vez el intervalo
en el
).
Distintas posibilidades que hay
que considerar.
Así, si representamos geométricamente las sucesivas aproximaciones, observamos que se obtiene una sucesión creciente acotada superiormente o una decreciente acotada inferiormente y, por tanto, en ambos casos la sucesión obtenida es convergente, como vamos a demostrar a continuación.
Demostración Es suficiente considerar una de las
cuatro posibilidades que hemos indicado anteriormente, por ejemplo la primera,
es decir, en la que y
,
. La función
es estrictamente creciente y
convexa en el intervalo considerado y, por lo tanto, si tomamos
, siendo
la única raíz de la ecuación
en
,
se tiene que
, con
; por
otro lado, al tratarse de una función convexa en dicho intervalo, se verifica
que
,
, donde
es la
recta tangente a la curva
en el punto
. En particular, se cumple
que
, deduciéndose el
único punto de corte de dicha recta
con el eje
, que es
.
Se obteniene de este modo una sucesión estrictamente decreciente (a no ser
que , para algún
)
y acotada inferiormente por
. De esta manera la
convergencia de la sucesión está asegurada hacia cierto valor
. Tomando límites en (metNR), se cumple que
, lo que equivale a
, ya
que
no puede anularse
por hipótesis. Concluimos pues que
, por la unicidad de la
raíz en
.
La ecuación (ecgen) puede escribirse de manera equivalente bajo multitud de expresiones que adopten una forma de tipo punto fijo
mediante la cual se puede generar
iterativamente una sucesión de valores a partir de un valor inicial
:
Geométricamente, se pasaría de la búsqueda de puntos de corte de la gráfica
con el eje de
abscisas
a la
búsqueda de puntos de intersección entre las dos gráficas de ecuaciones
e
. Pero precisamente, según
elijamos esta función
tendremos uno u otro método
de entre todos los posibles, siempre y cuando (ecptofijo) resulte equivalente a (ecgen). Unos serán convergentes y otros no.
Si es un punto fijo o raíz de la ecuación (ecptofijo), para el error de la
aproximación n-ésima
se tiene, en el supuesto de
que
sea derivable, que
para cierto
entre
y
.
Por tanto, el error en la etapa -ésima disminuye en valor
absoluto si la función
es menor que
en un entorno de la
solución y partimos de un punto de dicho entorno. Es decir, en esas condiciones
el método converge. Se tiene pues un resultado de convergencia local, ya que se están imponiendo restricciones respecto al
valor inicial
.
Así pues, cuando , existe un entorno centrado
en
en el cual hay
convergencia, al menos local, del método de iteración funcional (metiterfunc). La situación más favorable
se da cuando
. Y mejor
aún si algunas derivadas sucesivas siguen anulándose en
(digamos
, pues en tal caso
deduciéndose que
Se dice que el método es de orden
En la práctica, orden
(para
) significa que
cuando se está próximo a la solución en cada iteración el número de cifras
exactas se multiplica por
. De ahí el interés en que
el orden
sea mayor
que
. La dificultad en
tales métodos se reduce pues a encontrar una buena aproximación inicial. Para
ello, en una variable puede ser válido el uso del método de bisección; pero como
técnica general para la búsqueda de la aproximación inicial (válida tanto en una
como en varias variables) la conocida como método de
continuación es la más apropiada y se verá en el tema dedicado a los
sistemas de ecuaciones no lineales.
También podemos, basándonos en el concepto de contracción en espacios
métricos abstractos, dar un resultado de convergencia global para el método iterativo (metiterfunc), es decir partiendo de un
valor arbitrario de . El resultado teórico en
cuestión asegura la existencia de un único punto fijo para una aplicación contractiva
, definida en un espacio
métrico completo
con distancia
. Recordamos que se dice
que
es una contracción si se verifica la siguiente condición: existe una
constante
tal que
para cualesquiera
.
Así pues, basándonos en este teorema de punto fijo, podemos enunciar el
siguiente teorema de convergencia global.
Si es una función real
definida en cierto subconjunto cerrado (no necesariamente acotado)
de la recta
real, de manera que
,
,
y se tiene que
, con
, entonces
existe una única raíz de la ecuación
, que se puede obtener
como límite de la sucesión
obtenida mediante el
método iterativo (metiterfunc)
partiendo de un valor arbitrario
.
La demostración de este resultado es como sigue,
razonando por inducción. Se tiene que y, dado que
, vemos pues que
la sucesión
convergerá si y sólo si lo hace la serie
. Pero la desigualdad
(comparacion) y criterio de comparación de series nos
permiten escribir
que equivale a la convergencia
absoluta de la serie. Obtenemos como consecuencia la convergencia a
de la sucesión
. Por ser
cerrado, se cumple que
. Además,
(nótese que el hecho de ser
una aplicación
contractiva implica que, en particular, es continua).
La unicidad de la raíz se deduce también de la contractividad de , ya que, por reducción al
absurdo, si
y
fueran dos
puntos fijos de la misma, se tendría que
llegándose a una contradicción en
caso de que
.
Nótese que el hecho de que y que la función
fuese continua, ya
garantizaría la existencia de puntos fijos para esta función. No obstante, se
requeriría por ejemplo una condición de tipo Lipschitz como la impuesta en el
resultado anterior para poder asegurar además la unicidad de este punto fijo.
Esto se tendría asegurado por ejemplo si
con
,
. En este último caso se
tendría además que
y
Sea una función de clase
tal que la
ecuación
tenga
una solución
con,
con
. Sea
un valor
suficientemente próximo a la raíz. Entonces la sucesión construida a partir de
la iteración funcional (metiterfunc) convergerá hacia
con un orden de
convergencia
, cumpliéndose además que
La elección de la función para el método de Newton
Raphson es
siempre que exista
y no se
anule. Evidentemente, si para un valor
se anula la función
para dicho valor la función
vale precisamente
(
es punto fijo de
) y viceversa.
Puede probarse fácilmente que, si es dos veces derivable en
y
,
entonces
y, por tanto, la convergencia del
método de Newton-Raphson es al menos de orden dos (o cuadrática), lo que hace
que sea uno de los métodos más interesantes para resolver ecuaciones no
lineales. Cuando
es al menos de clase
podemos
calcular también la derivada segunda de
y, utilizando (ordenp), escribir
. Este límite puede ser
nulo en algunos casos (lo que implicaría una convergencia de orden superior).
El análisis de la convergencia de los métodos de la secante y de regula-falsi
no puede hacerse de esta misma forma, ya que ninguno de ellos puede expresarse
como un proceso de iteración funcional del tipo dado por (metiterfunc), al emplear ambos dos
iteraciones anteriores para el cálculo de la siguiente. No obstante, mediante
las correspondientes ecuaciones en diferencias (consultar
por ejemplo Kincaid-Cheney) se puede analizar el orden
de convergencia de ambos métodos. El de la secante resulta intermedio entre el
lineal y el cuadrático, , mientras que el de
regula-falsi no pasa de ser lineal (
).
Otra cuestión interesante a la hora de aplicar cualquiera de estos métodos
iterativos, dados por la relación (metiterfunc), sería la posibilidad de
acelerar la convergencia de los mismos, siempre que esto sea posible. Así, por
ejemplo, a partir de una cierta sucesión de valores aproximados se puede conseguir otra
sucesión
de la
siguiente manera, que constituye el llamado método de
Aitken:
Si la
sucesión
converge
a
más rápidamente que
, en el sentido de
que
(ver [KincaidCheney])
Vemos que en este método de aceleración de la convergencia intervienen tres
valores consecutivos de la sucesión; así pues, si partiéramos por ejemplo de los
tres primeros de la sucesión original,, y aplicásemos (metAitken) obteniendo
para calcular a continuación
y procediésemos
de esta manera reiteradamente, obtendríamos el conocido como método de Steffensen de aceleración de la convergencia; que no
debe confundirse con el también denominado método de
Steffensen de iteración funcional para la aproximación de raíces, que
consiste en tomar
y buscar sus puntos fijos. Se puede
demostrar que este método también converge cuadráticamente, como el de
Newton-Raphson.
D. Kincaid, W. Cheney, ``Análisis Numérico: Las matemáticas del cálculo científico'', Addison-Wesley Iberoamericana.
M. Gasca, ``Cálculo numérico: Resolución de ecuaciones y sistemas'', Librería Central (Zaragoza, 1987).
J.J. Quesada Molina, ``Ecuaciones Diferenciales, Análisis Numérico y Métodos Matemáticos'', Edit. Santa Rita (Monachil-Granada, 1996).