Resolución numérica de ecuaciones

Introducción

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 MATH para $f:$ MATH al menos continua y con distinto signo en los extremos de $I$.

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.

Métodos numéricos elementales

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, MATH. Lo interesante será que esta sucesión de aproximaciones converja hacia una raíz, pongamos $s\in\QTR{Bbb}{R}$, 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.

Método de bisección

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 $\left[ a,b\right] $. 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

Teorema

Sea MATH una función continua en $\left[ a,b\right] $ y tal que MATH MATH. Entonces existe MATH tal que MATH.

En la figura c02g1 se presenta la interpretación geométrica. Se considera el intervalo MATH. En los extremos la función cambia de signo, y hay más de un cero.
c02g1.wmf
Idea geométrica del teorema de Bolzano.

Descripción del algoritmo

Vemos pues cómo, si no se llega a encontrar la solución exacta $s=x_{k}$ para algún MATH, entonces se obtiene una sucesión de intervalos encajados MATH con MATH, de manera que, si MATH para todo MATH, entonces el error MATH tiende a cero cuando $k$ tiende a infinito, y, por la complitud del cuerpo de los números reales, el método resulta convergente, es decir MATH en $\QTR{Bbb}{R}$. 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 $tol$, en la aproximación; para ello debemos detener el proceso cuando la longitud del intervalo correspondiente sea menor o igual que $tol$; también podemos calcular de antemano el número de iteraciones resolviendo la inecuación MATH, 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 $a$ y $b$ para los extremos del intervalo considerado y $c$ para el punto medio del mismo. En la figura c02g2 se ilustra.
c02g2.wmf
Idea geométrica del método de bisección

ALGORITMO DE BISECCIÓN
Entrar $f$, $a$, $b$, $tol$
Mientras $b-a\geq tol$
Hacer MATH
Si MATH entonces $c$ es raíz. Fin
Si MATH entonces $b=c$.
Si MATH entonces $a=c$.

Método de regula-falsi

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 MATH; 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 MATH y MATH, $k\geq0$, y tomar la intersección de ésta con el eje $Ox$ (ver figura c02g3); llamemos también $c\equiv x_{k}$ 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, $[a_{k},c]$ o bien $[c,b_{k}]$. Bastará con repetir este proceso de forma recursiva para ir obteniendo cada vez mejores aproximaciones de la raíz buscada.
c02g3.wmf
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 $c$ como el punto de corte de la correspondiente recta secante, de ecuación MATH con el eje $Ox$: MATH 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 MATH de manera que MATH, siendo $s$ una raíz, que sabemos que posee la ecuación (ecgen) en el intervalo de partida.

Método de la secante

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 $x_{0}=a$ y $x_{1}=b$ mediante la siguiente fórmula (obtenida de la misma forma que en (metrfalsi), con las identificaciones MATHy $a_{k}=x_{k-1}$): MATH

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

Método de Whittaker

En el método de Whittaker, partiendo de un valor inicial $x_{0}$, se aproxima la raíz de la ecuación (ecgen) mediante la abscisa del correspondiente punto de corte con el eje $Ox$ de la recta que pasa por los sucesivos puntos MATH, $k\geq0$, y tiene como pendiente un valor prefijado MATH fijo (ver figura c02g4 ).
c02g4.wmf
Idea geométrica del método de Whittaker

Así pues, la ecuación de esas rectas es MATH y el cálculo de $x_{k+1}$ a partir de $x_{k}$ por MATH 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 $m$ que aproxime, para todas las iteraciones, el valor MATH de (metsec) y, por tanto, simplificar los cálculos de cada iteración.

Método de Müller

Este método se basa en la sustitución de la función $f$ en la ecuación (ecgen) por un polinomio de segundo grado que pase por tres puntos dados, MATH, MATH y MATH (con $k\geq2$), de la correspondiente gráfica MATH hallando posteriormente el punto de corte de dicha parábola con el eje $Ox$ mediante la conocida fórmula para este tipo ecuaciones de segundo grado (ver figura c02g5 ).
c02g5.wmf
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.

Método de Newton-Raphson

Deducción geométrica del algoritmo

En este caso se obtiene una sucesión de aproximaciones partiendo de un valor inicial $x_{0}$, que debe ser dado o elegido convenientemente. Una vez conocida una aproximación, digamos $x_{k}$, la siguien-te, $x_{k+1}$, se obtiene hallando el punto de corte de la correspondiente recta tangente a la curva de ecuación $y=f(x)$ en el punto MATH con el eje $Ox$ (ver figura c02g6 ).
c02g6.wmf
Idea geométrica del método de Newton-Raphson.

El proceso se repite sucesivamente, obteniéndose pues el siguiente método iterativo: para $x_{0}$ dado (bastará con tomar MATH en (metWhittaker), con la salvedad de que ahora esta pendiente irá cambiando en cada iteración) MATH

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 MATH, MATH es el método que acabamos de deducir aplicado a $f(x)=x^{2}-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 MATH

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 $x_{0}$ (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 $x_{0}$ del método de Newton-Raphson .

Hemos de hacer algunas observaciones al método. En primer lugar, $f$ 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 $x_{0}$ y $x_{1}$ para almacenar los valores de $x_{k}$ y $x_{k+1}$, respectivamente.

ALGORITMO DE NEWTON-RAPHSON
Entrar $f$, $x_{0}$, $tol$
Hacer MATH
Si MATH entonces escribir $x_{1}$. Fin
En caso contrario intercambiar $x_{0}$ por $x_{1}$ 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.

Proposición

Supongamos que la función $f$ admite derivada segunda en $\left[ a,b\right] $ y verifica las siguientes propiedades:

Entonces, partiendo de un valor inicial $x_{0}\in[a,b]$ para el que se

verifique que MATH, la sucesión de aproximaciones obtenida por el método de Newton-Raphson converge hacia la única raíz de la ecuación MATH en $\left[ a,b\right] $.

Antes de probar el teorema, observemos que bastará con tomar por ejemplo como $x_{0}$ el extremo del intervalo para el cual $f$ y $f^{\prime\prime}$ 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 $f$ 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 $Ox$ (mediante el cambio MATH) o al eje $Oy$ (mediante el cambio $x\rightarrow-x$, que transformaría a su vez el intervaloMATH en el MATH ).
c02g7.wmf
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 $f^{\prime}(x)>0$ y MATH, $\forall$ $x$ MATH. La función $f$ es estrictamente creciente y convexa en el intervalo considerado y, por lo tanto, si tomamos MATH, siendo $s$ la única raíz de la ecuación en $\left[ a,b\right] $, se tiene que MATH, con $f(x_{k})>0$; por otro lado, al tratarse de una función convexa en dicho intervalo, se verifica que MATH, $\forall$ $x$ MATH, donde MATH es la recta tangente a la curva $y=f(x)$ en el punto MATH. En particular, se cumple que MATH, deduciéndose el único punto de corte de dicha recta$\ t_{k}(x)$ con el eje $Ox$, que es MATH.

Se obteniene de este modo una sucesión estrictamente decreciente (a no ser que $x_{k+1}=s$, para algún $\ k\in\QTR{Bbb}{N}$) y acotada inferiormente por $s$. De esta manera la convergencia de la sucesión está asegurada hacia cierto valor MATH. Tomando límites en (metNR), se cumple que MATH, lo que equivale a $f(\widehat{s})=0$, ya que MATH no puede anularse por hipótesis. Concluimos pues que MATH, por la unicidad de la raíz en $\left[ a,b\right] $. $\Box$

Técnicas de iteración funcional

La ecuación (ecgen) puede escribirse de manera equivalente bajo multitud de expresiones que adopten una forma de tipo punto fijo

MATH mediante la cual se puede generar iterativamente una sucesión de valores a partir de un valor inicial $x_{0}$:

MATH

Geométricamente, se pasaría de la búsqueda de puntos de corte de la gráfica MATH con el eje de abscisas $Ox$ a la búsqueda de puntos de intersección entre las dos gráficas de ecuaciones $y=x$ e MATH. Pero precisamente, según elijamos esta función $g$ 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.

Orden de convergencia

Si $s\in\QTR{Bbb}{R}$ es un punto fijo o raíz de la ecuación (ecptofijo), para el error de la aproximación n-ésima $e_{k}=s-x_{k}$ se tiene, en el supuesto de que $g$ sea derivable, que

MATH para cierto $\xi$ entre $s$ y $x_{k}$.

Por tanto, el error en la etapa $k$-ésima disminuye en valor absoluto si la función $g^{\prime}$ es menor que $1$ 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 $x_{0}$.

Así pues, cuando MATH, existe un entorno centrado en $s$ 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 MATH. Y mejor aún si algunas derivadas sucesivas siguen anulándose en $s$ (digamos MATH, pues en tal caso

MATH deduciéndose que

MATH Se dice que el método es de orden $p\in\QTR{Bbb}{N}.$ En la práctica, orden $p$ (para $p>1$) significa que cuando se está próximo a la solución en cada iteración el número de cifras exactas se multiplica por $p$. De ahí el interés en que el orden $p$ sea mayor que $1$. 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 $x_{0}$. El resultado teórico en cuestión asegura la existencia de un único punto fijo para una aplicación contractiva MATH, definida en un espacio métrico completo MATH con distancia MATH. Recordamos que se dice que $T$ es una contracción si se verifica la siguiente condición: existe una constante MATH tal que MATH para cualesquiera $x,y\in\QTR{cal}{X}$ . Así pues, basándonos en este teorema de punto fijo, podemos enunciar el siguiente teorema de convergencia global.

Teorema

Si $g$ es una función real definida en cierto subconjunto cerrado (no necesariamente acotado) MATH de la recta real, de manera que MATH, $\,\forall x\in I$, y se tiene que MATH, con $L<1$, entonces existe una única raíz de la ecuación MATH, que se puede obtener como límite de la sucesión MATH obtenida mediante el método iterativo (metiterfunc) partiendo de un valor arbitrario $x_{0}\in I$.

Demostración

La demostración de este resultado es como sigue, razonando por inducción. Se tiene que MATH y, dado que MATH, vemos pues que la sucesión MATH convergerá si y sólo si lo hace la serie MATH. Pero la desigualdad (comparacion) y criterio de comparación de series nos permiten escribir MATH que equivale a la convergencia absoluta de la serie. Obtenemos como consecuencia la convergencia a $s$ de la sucesión MATH. Por ser $I$ cerrado, se cumple que $s\in I$. Además, MATH (nótese que el hecho de ser $g$ una aplicación contractiva implica que, en particular, es continua).

La unicidad de la raíz se deduce también de la contractividad de $g$, ya que, por reducción al absurdo, si$\ s_{1}$ y $s_{2}$ fueran dos puntos fijos de la misma, se tendría que MATH llegándose a una contradicción en caso de que $s_{2}\neq s_{1}$. $\Box$

Nótese que el hecho de que MATH y que la función $g$ 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 MATH con MATH, MATH. En este último caso se tendría además que MATH y MATH

Proposición

Sea $g$ una función de clase MATH tal que la ecuación MATH tenga una solución MATH con, MATH con MATH. Sea $x_{0}$ un valor suficientemente próximo a la raíz. Entonces la sucesión construida a partir de la iteración funcional (metiterfunc) convergerá hacia $s$ con un orden de convergencia $p=m$, cumpliéndose además que MATH

Deducción analítica del orden de convergencia del método de Newton-Raphson

La elección de la función $g$ para el método de Newton Raphson es

MATH siempre que exista $f^{\prime}$ y no se anule. Evidentemente, si para un valor $s$ se anula la función $f$ para dicho valor la función $g$ vale precisamente $s$ ($s$ es punto fijo de $g$) y viceversa.

Puede probarse fácilmente que, si $f$ es dos veces derivable en $s$ y $f^{\prime}(s)\neq0$, entonces MATH 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 $f$ es al menos de clase $C^{3}$ podemos calcular también la derivada segunda de $g$ y, utilizando (ordenp), escribir MATH. 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, MATH, mientras que el de regula-falsi no pasa de ser lineal ($p=1$).

Aceleración de la convergencia

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 MATH se puede conseguir otra sucesión MATH de la siguiente manera, que constituye el llamado método de Aitken:

MATH

Si MATH la sucesión MATH converge a $s$ más rápidamente que MATH, en el sentido de que MATH (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,MATH, y aplicásemos (metAitken) obteniendo MATHpara calcular a continuación MATH 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 MATHy buscar sus puntos fijos. Se puede demostrar que este método también converge cuadráticamente, como el de Newton-Raphson.

    1. D. Kincaid, W. Cheney, ``Análisis Numérico: Las matemáticas del cálculo científico'', Addison-Wesley Iberoamericana.

    2. M. Gasca, ``Cálculo numérico: Resolución de ecuaciones y sistemas'', Librería Central (Zaragoza, 1987).

    3. J.J. Quesada Molina, ``Ecuaciones Diferenciales, Análisis Numérico y Métodos Matemáticos'', Edit. Santa Rita (Monachil-Granada, 1996).