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 entonces es raíz. Fin |
Si entonces . |
Si entonces . |
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 entonces escribir . Fin |
En caso contrario intercambiar por y reiterar. |
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).