MARGARITA MATHEMATICA EN MEMORIA DE ´ JAVIER (CHICHO) GUADALUPE HERNANDEZ ´ JOSE (Luis Espa˜ nol y Juan L. Varona, editores), Servicio de Publicaciones, Universidad de La Rioja, Logro˜ no, Spain, 2001.

´ ´ ´ EL METODO DE HALLEY: POSIBLEMENTE, EL METODO MAS REDESCUBIERTO DEL MUNDO ´ A. EZQUERRO, JOSE ´ M. GUTIERREZ, ´ JOSE ´ MIGUEL A. HERNANDEZ Y M. AMPARO SALANOVA A nuestro amigo y compa˜ nero Chicho Abstract. Halley’s method is a third order iterative method to solve nonlinear equations. Many authors have deduced it independtly. Here we try to justify, by showing different constructions, why it is said that Halley’s method is the most “rediscovered” method in the mathematical literature.

Si tuvi´eramos que hacer una clasificaci´on con los m´as famosos m´etodos iterativos para resolver ecuaciones, no hay duda alguna de que el m´etodo de Newton ocupar´ıa el primer lugar. Seg´ un la opini´ on de varios autores y, teniendo en cuenta el n´ umero de publicaciones relacionadas con ´el, el m´etodo de Halley ocupar´ıa el segundo puesto en esta ficticia clasificaci´on. As´ı, por ejemplo, Traub, en su cl´ asica monograf´ıa [19, p´ ag. 91], dice: ✭✭El m´etodo de Halley puede compartir con el m´etodo de la secante la distinci´ on de ser el m´etodo m´as frecuentemente redescubierto en la literatura✮✮. Confirmando este hecho, Scavo y Thoo [17] presentan un excelente art´ıculo en el que se cuenta la historia del m´etodo, respaldada por numerosas referencias bibliogr´ aficas. Es costumbre extendida entre los investigadores el bautizar sus descubrimientos con su propio nombre o con el de un personaje relevante en la materia. La resoluci´ on num´erica de ecuaciones no lineales no escapa a esta costumbre. En este caso, dos de los m´etodos m´as conocidos, los m´etodos de Newton y Halley, reciben sus nombres de dos eminentes cient´ıficos brit´ anicos de finales del siglo XVII y comienzos del XVIII: Sir Isaac Newton y Edmund Halley. Pero, ¿hasta qu´e punto conoc´ıan estos autores los m´etodos cuya ✭✭paternidad✮✮ se les atribuye? Por lo que respecta al m´etodo de Newton, parece ser que la idea que pudo tener Sir Isaac dista bastante de lo que conocemos hoy en d´ıa. As´ı, en una carta a sus colegas Barrow y Collins en 1669, Newton mostraba un ejemplo para resolver num´ericamente una ecuaci´on. Newton ilustraba su t´ecnica con el ejemplo x3 − 2x − 5 = 0 y argumentaba de la siguiente manera: Por tanteo, se ve que la soluci´ on est´a cerca de 2. Haciendo x = 2 + ε y sustituyendo en la ecuaci´ on se obtiene: (1)

ε3 + 6ε2 + 10ε − 1 = 0.

2000 Mathematics Subject Classification. 65–02. Key words and phrases. Nonlinear equations, Halley’s method, third order methods. Este trabajo ha sido subvencionado por una ayuda de la DGES (ref. PB98-0198) y otra de la Universidad de La Rioja (ref. API-00/B16). 205

206

´ ´ J. A. EZQUERRO, J. M. GUTIERREZ, M. A. HERNANDEZ Y M. A. SALANOVA

Ignorando los t´erminos ε3 + 6ε2 con el pretexto de que ε es peque˜ no, se llega a que 10ε − 1 = 0 o´ ε = 0,1. Entonces x = 2,1 es una aproximaci´ on de la soluci´ on mejor que la inicial. Haciendo ahora ε = 0,1 + ν y sustituyendo en (1) se sigue que ν 3 + 6,3ν 2 + 11,23ν + 0,061 = 0. Ignorando de nuevo los t´erminos en ν de grado mayor o igual que dos, se llega a que ν = −0,0054 y, por tanto, x = 2,0946 es una aproximaci´on que mejora las anteriores. Newton indicaba que el proceso se puede repetir las veces que sean necesarias. Como vemos, la idea de Newton consiste en a˜ nadir un t´ermino corrector a una aproximaci´on inicial dada. Para obtener esta aproximaci´ on, lo que hace es truncar el binomio de Newton en el segundo t´ermino en expresiones del tipo (a + ε)n  an + nan−1 ε. De esta manera, para obtener el valor aproximado de ε simplemente hay que resolver una ecuaci´ on lineal. Escribiendo el problema con la notaci´ on actual y llamando p(x) = x3 − 2x − 5, tenemos que la nueva aproximaci´ on es 2−

1 p(2) =2+ = 2,1,  p (2) 10

que se corresponde con la conocida formulaci´ on del m´etodo de Newton (2)

xn+1 = xn −

p(xn ) . p (xn )

Ahora bien, no se tiene constancia de que Newton usara el c´ alculo diferencial ni de que expresara el proceso como un m´etodo iterativo en el sentido de que una aproximaci´on pueda ser considerada como punto de partida de la siguiente. La idea de iteraci´ on se atribuye a Joseph Raphson, quien adem´ as simplifica el aspecto operacional de la t´ecnica de Newton. En 1690 publica un tratado en el que se dan f´ ormulas expl´ıcitas para el t´ermino corrector para algunos casos particulares. En concreto, calcula los t´erminos correctores para las ecuaciones x3 − r = 0 y x3 − px − q = 0 que son, respectivamente, r − x30 3x20

y

q + px0 − x30 , 3x20 − p

on inicial. La contribuci´ on de Raphson ha sido tenida en siendo x0 la aproximaci´ cuenta hist´ oricamente, no en vano muchos autores denominan el proceso como m´etodo de Newton-Raphson. Sin embargo, en los trabajos de Raphson no se aprecia la conexi´on existente entre el t´ermino corrector, la funci´on que define la ecuaci´ on y su derivada. La incorporaci´ on del c´ alculo diferencial se debe a Thomas Simpson (1740) quien estableci´o el m´etodo tal y como lo conocemos actualmente, salvo aspectos notacionales (Simpson explicaba de forma ret´ orica c´omo obtener las aproximaciones sucesivas). Adem´as, Simpson extendi´ o el proceso a funciones cualesquiera, no solamente polinomios.

´ EL METODO DE HALLEY

207

Posteriormente, Lagrange (1798) expresa el m´etodo con la notaci´ on actual (2) y Fourier (1831) lo bautiza como ✭✭la m´ethode newtonienne✮✮. Quiz´ a, Fourier es el causante de la falta de reconocimiento para el trabajo de Simpson. Con estas l´ıneas hemos querido indicar que, en ocasiones, la ✭✭paternidad✮✮ de un m´etodo no se debe exclusivamente al cient´ıfico que le da el nombre, sino que es contribuci´ on de varios investigadores. Para profundizar m´ as en la historia del m´etodo de Newton, recomendamos el trabajo de Ypma [20]. Algo parecido ocurre con el m´etodo de Halley. Siguiendo a Bailey [1], podemos ver que Halley qued´ o impresionado por los trabajos del matem´ atico franc´es Thomas Fautet de Lagny, quien en 1691 dio f´ormulas (racionales e irracionales) para aproximar la ra´ız c´ ubica de un n´ umero. Halley, en 1708 escribe en la revista londinense Miscellanea Curiosa lo siguiente: ✭✭Habiendo comprobado la eficiencia de estas f´ ormulas (las de de Lagny), voy a intentar encontrar la demostraci´ on✮✮. Y he aqu´ı la demostraci´on de Halley para la f´ ormula racional. Se trata de aproximar la ra´ız c´ ubica de un n´ umero natural escrito en la forma a3 + b, con a, b ∈ Z. Para ello, Halley busca un t´ermino corrector ε de forma que (a + ε)3 = a3 + b. De aqu´ı se deduce que b = 3a2 ε + 3aε2 + ε3 . Despreciando el t´ermino ε3 , se llega a que b  3a2 ε + 3aε2 y se obtienen las siguientes aproximaciones para ε: ε2 b ε+  ε, 2 3a a

(3) (4)

3a2

b  ε. + 3aε

Sustituyendo (3) en (4) se llega a b

ab = 3 . b 3a +b 3a2 + 3a 2 3a Esta es la f´ ormula racional dada por de Lagny:
ab 3 a3 + b  a + 3 . 3a + b Si escribimos con notaci´on actual el problema, llamando f(x) = x3 − a3 − b y tomando x0 = a, la expresi´ on anterior se corresponde efectivamente con el primer paso del m´etodo de Halley: 2 f(xn ) (5) xn+1 = xn − , 2 − Lf (xn ) f  (xn ) ε

donde f(xn )f  (xn ) . f  (xn )2 Parece ser que Halley intent´o extender sus m´etodos para obtener aproximaciones incluso para ecuaciones dadas por funciones trascendentes. De hecho, en un trabajo suyo de 1694 aparecen impl´ıcitamente aproximaciones del tipo (5) y tambi´en aproximaciones dadas por f´ormulas irracionales. Lo que no est´ a claro es d´onde aparecieron por primera vez estas f´ ormulas expl´ıcitamente. As´ı, debemos esperar un siglo (6)

Lf (xn ) =

208

´ ´ J. A. EZQUERRO, J. M. GUTIERREZ, M. A. HERNANDEZ Y M. A. SALANOVA

y medio hasta que en un trabajo de Schr¨ oder (1870) podemos encontrar la iteraci´ on de Halley tal y como la conocemos actualmente. No obstante, en dicho trabajo, el m´etodo de Halley s´olo aparece de pasada y sin ninguna referencia al propio Halley. Posteriormente, el m´etodo de Halley ir´ a apareciendo con cierta frecuencia en los trabajos de muchos autores que lo ir´ an ✭✭redescubriendo✮✮. Para hacernos una idea de este hecho, podemos consultar la historia del m´etodo que aparece en el trabajo de Scavo y Thoo [17]. En lo que sigue, vamos a exponer diferentes formas de ✭✭descubrir✮✮ el m´etodo de Halley. Pretendemos exponer estas t´ecnicas por el inter´es que tienen por s´ı mismas y porque, adem´ as, estos procedimientos nos permitir´an, mediante ligeras variantes, construir otros procesos iterativos. No se pretende asignar un descubridor a cada una de estas t´ecnicas. En consecuencia, entendemos que los autores que citamos han trabajado en ellas, pero, eso s´ı, no tienen porqu´e ser los primeros ni los u ´ nicos en hacerlo. No es nuestra intenci´on realizar aqu´ı un estudio de las caracter´ısticas de los m´etodos que tratamos (convergencia local o global, orden de convergencia, an´ alisis del error, etc.). Simplemente queremos mostrar distintas maneras de construir procesos iterativos y, en particular, el m´etodo de Halley. Del mismo modo, tampoco precisaremos las condiciones que han de cumplirse para que estas construcciones puedan realizarse. As´ı, supondremos que las funciones que manejaremos son todo lo continuas y derivables que deban ser. Primera forma: a partir del desarrollo de Taylor. Seg´ un se puede leer en el art´ıculo de Bateman [2], el trabajo de Halley para extender su m´etodo para aproximar soluciones a ecuaciones trascendentes inspir´o a Brook Taylor en la consecuci´ on de su afamado teorema. El propio Taylor escribe en una carta a un colega en 1712: ✭✭Estaba pensando en c´ omo aplicar el m´etodo del Dr. Halley a la ecuaci´ on de Kepler cuando se me ocurri´o la forma de extender su m´etodo de extracci´ on de ra´ıces a todos los problemas . . . ✮✮ La idea que presentamos a continuaci´on es, probablemente, la manera m´ as habitual de ✭✭redescubrir✮✮ el m´etodo de Halley. En ella empleamos el desarrollo de Taylor de una funci´ on, que como acabamos de ver es un resultado que se public´o con posterioridad a los trabajos de Newton, Raphson o Halley. Partiendo de una aproximaci´ on x0 de la soluci´ on de una ecuaci´ on f(x) = 0, buscamos una aproximaci´ on mejor de la forma x0 +ε. Lo ideal ser´ıa que f(x0 +ε) = 0. Teniendo en cuenta esto y usando el desarrollo de Taylor, llegamos a (7)

0 = f(x0 + ε) = f(x0 ) + f  (x0 )ε +

f  (x0 ) 2 f  (x0 ) 3 ε + ε +··· 2 6

La idea de Newton o Halley de truncar el binomio de Newton se extiende ahora al desarrollo de Taylor. De esta manera, truncando (7) en el segundo sumando se llega al m´etodo de Newton. Si a˜ nadimos un t´ermino m´as, es decir, si truncamos (7) en el t´ermino con potencias en ε2 , se obtiene (8)

0  f(x0 ) + f  (x0 )ε +

f  (x0 ) 2 ε , 2

´ EL METODO DE HALLEY

209

y de aqu´ı, f(x) 2
. 1 ± 1 − 2Lf (x) f  (x) Para evitar restar dos cantidades similares, se recomienda [19] elegir la ra´ız cuadrada con signo positivo. Por tanto, la nueva aproximaci´ on pasa a ser f(x0 ) 2
(9) x1 = x0 − . 1 + 1 − 2Lf (x0 ) f  (x0 ) ε

Este m´etodo, conocido como m´etodo de Euler o m´etodo de Halley irracional ha sido estudiado, entre otros, por Traub [19, p´ ag. 93] y Melman [13]. La expresi´on correspondiente para el c´alculo de ra´ıces c´ ubicas, es decir, cuando f(x) = x3 − a3 − b era conocida por Halley (y por de Lagny):  a2 a b + , (x0 = a). x1 = + 2 2 3a Adem´as parece ser que Halley era m´as partidario de las expresiones que depend´ıan del c´alculo de ra´ıces cuadradas que de las que conten´ıan fracciones con denominadores muy grandes. Con este comentario, Halley anuncia con siglos de antelaci´ on uno de los problemas con los que nos podemos encontrar al trabajar num´ericamente. Halley cre´ıa incluso que la f´ ormula irracional era m´ as precisa que la racional1. Sin embargo, no es este el m´etodo que ha heredado su nombre y que, con el transcurso del tiempo, ha adquirido m´ as fama. El m´etodo en cuesti´on evita el c´alculo de ra´ıces cuadradas. Para ello, en la aproximaci´ on (8) se hace f(x0 ) . f  (x0 ) f  (x0 ) + ε 2 Para no tener que resolver una ecuaci´ on de segundo grado y no tener que calcular ra´ıces cuadradas, se sustituye el t´ermino corrector ε que aparece en el denominador anterior por la correcci´ on del m´etodo de Newton, ε = −f(x0 )/f  (x0 ), para obtener as´ı la correcci´on de Halley: f(x0 ) . ε− f(x0 )f  (x0 ) f  (x0 ) − 2f  (x0 ) (10)

ε−

El m´etodo resultante es el que se conoce como m´etodo de Halley (5): xn+1 = xn −

2 f(xn ) . 2 − Lf (xn ) f  (xn )

Se observa el paralelismo que hay entre la forma de deducir el m´etodo que acabamos de comentar y el trabajo original de Halley para el c´ alculo de ra´ıces c´ ubicas explicado con anterioridad. Extensiones de esta misma t´ecnica han permitido encontrar gran n´ umero de procesos iterativos. Si en lugar de sustituir el t´ermino corrector ε que aparece en el 1Esta afirmaci´ on no es cierta en general, como lo demuestra el contraejemplo dado por f (x) = ex − 2, con x0 = 1,3.

210

´ ´ J. A. EZQUERRO, J. M. GUTIERREZ, M. A. HERNANDEZ Y M. A. SALANOVA

denominador de (10) por la correcci´ on del m´etodo de Newton se sustituye por otras correcciones, se obtienen otros procesos. Esta t´ecnica, denominada de reemplazamiento, ha sido utilizada por autores como Neta [14] para definir una treintena de nuevos procesos iterativos de tercer orden. Stewart [18] y Hamilton [9], en lugar de trabajar con la ecuaci´ on (7), la transforman en un sistema de dos ecuaciones como sigue:  f  (x0 ) 2 f  (x0 ) 3   f  (x0 )ε + ε = −f(x0 ) − ε +··· 2 6    f(x )ε + f  (x )ε2 = − f (x0 ) ε3 + · · · 0 0 2 k Despreciando los t´erminos en ε con k ≥ 3, denotando u y v a las aproximaciones de ε y ε2 respectivamente, y resolviendo el sistema resultante por la regla de Cramer, se tiene que    −f(x0 ) f  (x0 )/2     0 f  (x0 )  2 f(x0 )  = ε  u =   ,   (x )  2 − L (x ) f f (x ) f (x )/2 f 0 0 0 0    f(x0 ) f  (x0 )  que es la correcci´on del m´etodo de Halley (5). Extensiones de esta t´ecnica permiten encontrar procesos iterativos de distintos o´rdenes, a partir de expresiones que dependen de determinantes. Para obtener una informaci´ on m´ as detallada, se puede consultar el trabajo de Kalantari, Kalantari y Zaare-Nahandi [11]. Segunda forma: construcci´ on geom´ etrica. Es bien conocido que el m´etodo de Newton tiene una clara interpretaci´ on geom´etrica, ya que, dado una aproximaci´ on inicial x0 , la siguiente aproximaci´ on, x1 , se obtiene como la intersecci´on de la recta tangente a la funci´ on f(x) en el punto x0 con el eje de abscisas. As´ı, x1 es la intersecci´on de las rectas y − f(x0 ) = f  (x0 )(x − x0 ) e y = 0, es decir, f(x0 ) . f  (x0 ) Por este motivo, al m´etodo de Newton se le suele denominar tambi´en como m´etodo de la tangente. El m´etodo de reemplazamiento introducido en la secci´ on anterior no nos permite dar una interpretaci´ on geom´etrica clara sobre la construcci´ on del m´etodo de Halley. No obstante, en 1952 Salehov [16] sugiere que el m´etodo de Halley puede obtenerse aproximando una funci´ on por medio de hip´erbolas. Siguiendo a Scavo y Thoo [17], para aproximar la soluci´ on de una ecuaci´ on f(x) = 0, a partir de un valor inicial x0 , vamos a buscar una hip´erbola de la forma (x − x0 ) + c y(x) = a(x − x0 ) + b x1 = x0 −

que cumpla y(x0 ) = f(x0 ), y (x0 ) = f  (x0 ) e y (x0 ) = f  (x0 ). Este tipo de curvas, denominadas osculatrices, coinciden con f en x0 hasta la segunda derivada. En consecuencia, estas hip´erbolas aproximan mejor a una curva que la recta tangente, donde la coincidencia s´ olo es hasta la primera derivada. Esta mejor aproximaci´ on se

´ EL METODO DE HALLEY

211

ve reflejada en una convergencia m´ as r´apida para el m´etodo de Halley que para el de Newton. El primero tiene orden tres mientras que el segundo tiene orden dos. De forma breve, se puede decir que el orden de convergencia mide la velocidad con la que una sucesi´on se aproxima a su l´ımite: a mayor orden, mayor velocidad. Volveremos a tratar el tema del orden en la u ´ltima secci´on. Para que la hip´erbola cumpla las condiciones requeridas sus coeficientes a, b y c deben satisfacer: c  = f(x0 ),   b    b − ac = f  (x0 ), 2 b       2a(ac − b) = f  (x0 ), b3 es decir,  −f(x0 )   , a=   2  2f (x0 ) − f(x0 )f  (x0 )     2f  (x0 ) , b=   2f (x0 )2 − f(x0 )f  (x0 )     2f(x0 )f  (x0 )   . c=  2f (x0 )2 − f(x0 )f  (x0 ) Una vez encontrada la hip´erbola osculatriz, tomamos como nueva aproximaci´on su intersecci´on con el eje de abscisas, es decir el punto tal que y(x1 ) = 0. Dada la forma de la hip´erbola, se tiene que 2f(x0 )f  (x0 ) x1 = x0 − c = x0 −  2f (x0 )2 − f(x0 )f  (x0 ) 2 f(x0 ) = x0 − . 2 − Lf (x0 ) f  (x0 ) Debido a esta interpretaci´ on geom´etrica, el m´etodo de Halley es llamado a veces el m´etodo de las hip´erbolas tangentes. Ideas similares se pueden utilizar para construir otros procesos iterativos, aproximando una funci´ on por otro tipo de curvas osculatrices. Por ejemplo, si en lugar de hip´erbolas consideramos par´ abolas de la forma y(x) = a(x − x0 )2 + b(x − x0 ) + c, se obtiene el m´etodo de Euler (9). Si las par´ abolas son de la forma y(x)2 + ay(x) + b(x − x0 ) + c = 0, el m´etodo obtenido es el de Chebyshev, m´etodo que presentamos en la siguiente secci´on. Tercera forma: a partir de la interpolaci´ on racional inversa. Presentamos en esta secci´on otra t´ecnica para construir procesos iterativos, conocida como interpolaci´ on inversa. En breve, podemos resumir el problema diciendo que se trata de resolver una ecuaci´on f(x) = 0. Para ello, consideramos la funci´ on que a cada x le hace corresponder y = f(x), y denotamos por x = φ(y) a su funci´ on inversa.

212

´ ´ J. A. EZQUERRO, J. M. GUTIERREZ, M. A. HERNANDEZ Y M. A. SALANOVA

El problema de aproximar la soluci´ on de f(x) = 0 se transforma entonces en el de aproximar φ(0). Se nos ofrecen ahora distintas posibilidades, en funci´ on de las diferentes maneras de aproximar la funci´ on φ. La primera de ellas consiste en aproximar φ(0) usando el desarrollo de Taylor hasta el orden dos. Esta idea se atribuye al matem´ atico ruso Pafnuty Lvovich Chebyshev, quien present´ o el m´etodo que lleva su nombre en un concurso estudiantil celebrado entre 1840 y 1841 y en el que consigui´ o la medalla de plata. Aunque existe una cierta unanimidad entre los cient´ıficos en atribuir a Chebyshev la ✭✭paternidad✮✮ del m´etodo, hay quien se lo atribuye a Leonhard Euler, lo cual no ser´ıa de extra˜ nar, debido a la prolijidad cient´ıfica de este insigne matem´atico y a que Chebyshev conoc´ıa bien la obra de Euler. No en vano, Chebyshev contribuy´ o a preparar una edici´ on con algunas de las obras de Euler a finales de 1840. Volviendo al problema de la interpolaci´ on inversa, supongamos que x0 es una aproximaci´on inicial de la soluci´ on de f(x) = 0 y sea y0 = f(x0 ). Entonces 1 φ(0)  φ(y0 ) − φ (y0 )y0 + φ (y0 )y02 2 f(x0 ) 1 f(x0 )2 f  (x0 ) = x0 −  − f (x0 ) 2 f  (x0 )3  1 f(x0 ) = x0 − 1 + Lf (x0 ) , 2 f  (x0 ) donde Lf (x) se define en (6). A la iteraci´on resultante  1 f(xn ) (11) xn+1 = xn − 1 + Lf (xn ) 2 f  (xn ) se la conoce como m´etodo de Chebyshev, m´etodo de la interpolaci´ on cuadr´ atica inversa o m´etodo de las par´ abolas tangentes. Desde otro punto de vista, podemos expresar la idea de Chebyshev indicando que se ha aproximado la funci´ on inversa por un polinomio de segundo grado φ(y)  a + b(y − y0 ) + c(y − y0 )2 = p(y), eligiendo los coeficientes de forma que φ(k)(y0 ) = p(k)(y0 ) para k = 0, 1, 2. Evidentemente, este polinomio es el polinomio de Taylor de segundo orden. Pero se pueden buscar otro tipo de funciones de aproximaci´ on, como por ejemplo, aproximantes de Pad´e, que es lo que hace Traub en [19, Cap. 5]. Si aproximamos φ(y) por una funci´ on racional del tipo φ(y) 

(y − y0 ) + a = R(y), b(y − y0 ) + c

exigiendo que φ(k)(y0 ) = R(k)(y0 ) para k = 0, 1, 2, se llega a que R(y) =

(2f  (x0 ) + x0 f  (x0 ))(y − y0 ) + 2x0 f  (x0 )2 . f  (x0 )(y − y0 ) + 2f  (x0 )2

´ EL METODO DE HALLEY

213

Eligiendo como nueva aproximaci´ on el valor x1 = R(0) se obtiene x1 =

−f(x0 )(2f  (x0 ) + x0 f  (x0 )) + 2x0 f  (x0 )2 , −f(x0 )f  (x0 ) + 2f  (x0 )2

expresi´on que, despu´es de operar adecuadamente, se puede escribir en la conocida forma 2 f(x0 ) x1 = x0 − , 2 − Lf (x0 ) f  (x0 ) es decir, de nuevo, el m´etodo de Halley (5). Cuarta forma: usando fracciones continuas. El procedimiento que mostramos en esta secci´on es una variante de la interpolaci´ on inversa introducida en la secci´ on precedente. La idea fue presentada por Frame en [5]. Manteniendo la misma notaci´ on de la secci´on anterior, es decir, denotando x = φ(y) a la funci´ on inversa de y = f(x), Frame obtuvo el siguiente desarrollo en serie de potencias para el t´ermino corrector ε:

ε = r − x0 = φ(0) − φ(y0 ) = ck ν k k≥1

f  (x0 ) 2 3f  (x0 ) − f  (x0 )f  (x0 ) 3 =ν−  ν + ν +··· , 2f (x0 ) 6f  (x0 ) donde r denota a la soluci´ on de f(x) = 0, ν = −f(x0 )/f  (x0 ) y c1 = 1, c2 =   −f (x0 )/(2f (x0 )), etc. En la secci´ on anterior hemos aproximado la funci´ on inversa por polinomios o aproximantes de Pad´e. En esta secci´on vamos a desarrollar en fracciones continuas el desarrollo en serie de potencias del t´ermino corrector:

a1 (12) ε= ck ν k = b 0 + a2 b1 + k≥1 a3 b2 + b3 + · · · aj bj−1 + bj + ρj+1 aj donde ρj = . bj + ρj+1 Dependiendo de las diversas fracciones continuas desarrolladas y del lugar en el que las trunquemos, obtendremos diversos procesos iterativos. En el desarrollo propuesto por Frame [5] se tiene que b0 = 0, bj = 1 ∀j ≥ 1, a1 = ν y aj+1 es el primer t´ermino del desarrollo en serie de potencias de ν del cociente aj − ρj ρj+1 = , ρj  donde ρ1 = ε = k≥1 ck ν k . As´ı, si truncamos la fracci´ on continua en el primer t´ermino, tenemos que ε  a1 = ν = −

f(x0 ) f  (x0 )

y la correcci´on que se obtiene es la del m´etodo de Newton.

214

´ ´ J. A. EZQUERRO, J. M. GUTIERREZ, M. A. HERNANDEZ Y M. A. SALANOVA

Para truncar en el segundo t´ermino se necesita calcular a2 que, como hemos dicho, es el primer t´ermino del desarrollo en serie de potencias de ρ2 . Como a1 − ρ1 ρ2 = = −c2 ν + (c22 − c3 )ν 3 + · · · , ρ1 se tiene que a2 = −c2 ν = −Lf (x0 )/2. Por tanto, ε

a1 2 f(x0 ) =− , 1 + a2 2 − Lf (x0 ) f  (x0 )

que es la correcci´on del m´etodo de Halley (5). Si truncamos la fracci´ on continua en t´erminos superiores, se obtienen m´etodos de orden superior, aunque el propio Frame advierte que, siguiendo la t´ecnica indicada unas l´ıneas m´as arriba, los c´ omputos se complican excesivamente a partir del cuarto t´ermino. No obstante, en el mismo art´ıculo [5] se proporciona otro procedimiento para el c´ alculo de los numeradores aj y los restos ρj mediante el uso de determinantes. Si en lugar del desarrollo en fracciones continuas considerado por Frame, se tienen en cuenta otros desarrollos, se obtienen otros procesos iterativos. Por ejemplo, siguiendo a Kincaid y Cheney ag. 417], se puede convertir una serie en fracci´ on [12, p´ continua. De hecho, la serie k≥1 ck ν k admite el desarrollo (12) con b0 = 0, a1 = 1, b1 = 1/ν, 2  1 1 1 aj = − , bj = + , j ≥ 2. cj−1 ν j−1 cj−1 ν j−1 cj ν j Si truncamos la fracci´ on continua en el segundo t´ermino, el t´ermino corrector que se obtiene es el correspondiente al m´etodo de Chebyshev (11). Quinta forma: como aceleraci´ on del m´ etodo de Newton. El problema que intentan resolver los procesos iterativos es el de aproximar la soluci´on de una ecuaci´ on f(x) = 0. En ocasiones, puede resultar m´ as sencillo resolver una ecuaci´on modificada f(x)g(x) = 0, donde g(x) es una funci´ on convenientemente elegida. Entre otras cosas, deber´a cumplir que no se anule en un entorno de la ra´ız buscada: si r es tal que f(r) = 0, entonces g(x) = 0 en un entorno de r. Si para aproximar la soluci´ on de una ecuaci´ on empleamos el m´etodo de Newton, la situaci´ on ideal ser´ıa que la funci´ on que manejamos sea una recta, ya que, evidentemente, en este caso la soluci´on exacta se alcanzar´ıa en la primera iteraci´ on. Pero bueno, en la pr´ actica no tiene mucho sentido usar el m´etodo de Newton para resolver una ecuaci´on lineal. De hecho, los m´etodos iterativos se usan para resolver ecuaciones no lineales. Pero a partir de esta idea y como sugiere Brown en [3], podemos usar la funci´ on g(x) para ✭✭linealizar✮✮ lo m´as posible la ecuaci´on f(x)g(x) = 0. ¿En qu´e sentido?; por ejemplo, exigiendo que (f(x)g(x)) se anule en x = r. Este ✭✭aumento de la linealidad✮✮ de la ecuaci´on en un entorno de r deber´ıa repercutir en la velocidad de convergencia del m´etodo de Newton. Recogiendo esta idea, Gerlach [7] prueba que cuando se aplica el m´etodo de Newton a una ecuaci´ on Ψ(x) = 0 donde la soluci´ on r cumple que Ψ (r) = 0, Ψ (r) =  0 y Ψ (r) = 0, entonces se alcanza la convergencia de orden tres, cuando en general,

´ EL METODO DE HALLEY

215

el orden de convergencia del m´etodo de Newton es dos. Al tener mayor orden, la sucesi´on obtenida para las funciones Ψ converge m´ as r´apidamente. En la u ´ltima secci´on trataremos con un poco m´ as de detalle el tema del orden de convergencia. En resumen, si Ψ(x) = f(x)g(x), buscamos g de manera que Ψ (r) = 0, es decir, f  (r)g(r) + 2f  (r)g (r) = 0. Despu´es de resolver esta ecuaci´on diferencial se llega a que la soluci´ on es (salvo constantes multiplicativas) 1 g(x) =
. f  (x) En consecuencia la idea de Brown
y Gerlach consiste en aplicar el m´etodo de Newton a la funci´ on Ψ(x) = f(x)/ f  (x) en lugar de a f(x). ¿Cu´ al es la expresi´on del m´etodo resultante? Ψ(x) 2 f(xn ) xn+1 = xn −  = xn − , Ψ (x) 2 − Lf (xn ) f  (xn ) expresi´on que, como no pod´ıa ser menos, coincide con la del m´etodo de Halley (5). Sexta forma: a partir de iteraciones de Laguerre. El m´etodo de Laguerre es un conocido m´etodo para encontrar las ra´ıces de un polinomio. De forma resumida, expondremos aqu´ı en qu´e consiste para mostrar despu´es una generalizaci´on que nos llevar´ a, entre otros, al m´etodo de Halley. Para un estudio m´ as detallado del m´etodo de Laguerre, recomendamos al lector el libro de Ostrowski [15, Ap´endice O]. Para ver c´ omo se construye la iteraci´on de Laguerre, vamos a considerar un polinomio p(x) de grado n con ra´ıces reales x1 ≤ x2 ≤ · · · ≤ xn : p(x) =

n (x − xi ). i=1

La idea del m´etodo de Laguerre parte de la derivaci´ on logar´ıtmica de p(x). As´ı, denotando ai = 1/(x − xi ), A = p (x)/p(x) y B = −dA/dx se llega a A = (log p(x)) =

n

i=1

1 ai , = (x − xi ) i=1 n

d(p (x)/p(x)) 1 p (x)2 − p(x)p (x)

B =− = = a2i . = 2 2 dx p(x) (x − xi ) n

n

i=1

i=1

La desigualdad de Cauchy-Schwarz nos permite deducir la siguiente desigualdad: (A − a1 )2 = (a2 + · · · + an )2 ≤ (n − 1)(a22 + · · · + a2n ) = (n − 1)(B − a21 ), o equivalentemente na21 − 2Aa1 − (n − 1)B + A2 ≤ 0. Entonces, a1 se encuentra situado entre las ra´ıces de la ecuaci´on (13)

nu2 − 2Au − (n − 1)B + A2 = 0,

216

´ ´ J. A. EZQUERRO, J. M. GUTIERREZ, M. A. HERNANDEZ Y M. A. SALANOVA

que son



(n − 1)(nB − A2 ) A + (n − 1)(nB − A2 ) , u2 = . u1 = n n Es m´as, se tiene que a2 , . . . , an est´an entre u1 y u2 , es decir, A−

1 ≤ u2 , x − xj

u 1 ≤ aj =

(14)

j = 1, . . . , n.

Este hecho permite encontrar mejores aproximaciones de las ra´ıces de un polinomio a partir de una aproximaci´ on dada. En efecto, si la aproximaci´on inicial x est´a entre dos ra´ıces consecutivas, digamos xj < x < xj+1 , entonces

1 1 1 si |F (x) − r| = O(|x − r|p).

218

´ ´ J. A. EZQUERRO, J. M. GUTIERREZ, M. A. HERNANDEZ Y M. A. SALANOVA

Es sencillo comprobar que esta condici´on de orden se satisface si F (x) posee en un entorno de r derivadas hasta el orden p y, adem´ as, F (r) = r, F  (r) = F  (r) = · · · = F (p−1)(r) = 0,

F (p)(r) = 0.

Entonces el proceso iterativo dado por la funci´ on de iteraci´ on F tiene exactamente orden p. Esto quiere decir que en cada iteraci´on el numero de cifras significativas con las que se va aproximando la soluci´ on r se multiplica aproximadamente por p. Tanto la definici´ on de orden como la caracterizaci´on que hemos dado se atribuyen al matem´atico de finales del siglo XIX E. Schr¨ oder. Para profundizar un poco m´ as en este concepto, se puede consultar el art´ıculo de Ehrmann [4]. A partir de la condici´ on suficiente anterior, y a la vista de la forma que tienen algunos de los procesos iterativos m´as conocidos, Gander [6] se plantea el estudio de los m´etodos dados por la siguiente funci´ on de iteraci´ on: (18)

F (x) = x −

f(x) G(x). f  (x)

En concreto, el objetivo es caracterizar la funci´ on G para que el proceso definido por la funci´ on anterior tenga orden tres. Denotamos u(x) = f(x)/f  (x). Notemos que u (x) = 1 − Lf (x), con Lf (x) definido en (6). Como F (r) = r, para que (18) defina un proceso iterativo de segundo orden, bastar´ a con que F  (r) = 1 − u (r)G(r) − u(r)G(r) = 1 − G(r) = 0, es decir, G(r) = 1. Como caso particular, para la funci´ on constante G(x) = 1, obtenemos el m´etodo de Newton. Si queremos m´etodos con convergencia c´ ubica, necesitamos a˜ nadir la condici´ on F (r) = 0, que se traduce en f  (r) G(r) − 2G (r) = 0. f  (r) Teniendo en cuenta que G(r) = 1, nos queda G (r) =

f  (r) . 2f  (r)

Estas condiciones no nos ayudan demasiado a elegir apropiadamente una funci´ on G, puesto que dependen de la ra´ız r. Sin embargo, si consideramos G(x) = H(Lf (x)), estas condiciones se traducen en G(r) = H(0) = 1, G (r) = H  (0)Lf (r) = H  (0)

f  (r) f  (r) = . f  (r) 2f  (r)

Esto nos permite construir f´ acilmente funciones H que cumplan 1 H(0) = 1, H  (0) = . 2

´ EL METODO DE HALLEY

219

Algunos ejemplos de estas funciones nos han ido apareciendo a lo largo de este trabajo: 2

El m´etodo de Halley (5): H(t) = (1 − 2t )−1 = 1 + 2t + t4 + · · · El m´etodo de Chebyshev (11): H(t) = 1 + 2t √ 2 El m´etodo de Euler (9): H(t) = 2(1 + 1 − 2t)−1 = 1 + 2t + t2 + · · · 2 El m´etodo de Ostrowski (17): H(t) = 2(1 − t)−1/2 = 1 + 2t + 3t8 + · · · La familia de m´etodos de Hansen y Patrick (16): −1


t (α + 3)t2 =1+ + +··· H(t) = (α + 1) α + 1 − (α + 1)t) 2 8 Otra familia de m´etodos de tercer orden que tambi´en se ajusta a la forma (18) es la estudiada en [8]. En este caso, 1. 2. 3. 4. 5.

t λt2 t =1+ + +··· 2(1 − λt) 2 2 siendo λ un par´ ametro real. Esta familia incluye, entre otros, un m´etodo muy conocido y que no nos hab´ıa aparecido hasta ahora. Este es el denominado m´etodo super-Halley y se corresponde con el valor λ = 1. Como conclusi´on, podemos decir que muchos de los procesos iterativos de tercer orden se ajustan a la forma de la funci´ on de iteraci´ on (18). Muchos, pero no todos. Se puede probar [4] que todos los procesos iterativos de tercer orden est´ an dados por la funci´ on de iteraci´ on: H(t) = 1 +

G(x) = H(Lf (x)) + f(x)2 b(x), donde b es una funci´ on arbitraria que est´ a acotada para valores pr´ oximos de la ra´ız r. Referencias [1] D. F. Bailey, A historical survey of solution by functional iteration, Mathematics Magazine 62 (1989), 155–166. [2] H. Bateman, Halley’s method for solving equations, Amer. Math. Monthly 45 (1938), 11–17. [3] G. H. Brown Jr., On Halley’s variation of Newton’s method, Amer. Math. Monthly 84 (1977), 726–728. [4] H. Ehrmann, Konstruktion und Durchf¨ uhrung von Iterationsverfahren h¨ oherer Ordn¨ ung, Arch. Rational Mech. Anal. 4 (1959), 65–88. [5] J. S. Frame, The solution of equations by continued fractions, Amer. Math. Monthly 60 (1953), 293–305. [6] W. Gander, On Halley’s iteration method, Amer. Math. Monthly 92 (1985), 131–134. [7] J. Gerlach, Accelerated convergence in Newton’s method, SIAM Rev. 36 (1994), 272–276. [8] J. M. Guti´errez y M. A. Hern´ andez, A family of Chebyshev-Halley type methods in Banach spaces, Bull. Austral. Math. Soc. 55 (1997), 113–130. [9] H. J. Hamilton, A type of variation of Newton’s method, Amer. Math. Monthly 57 (1950), 517–522. [10] E. Hansen y M. Patrick, A family of root finding methods, Numer. Math. 27 (1977), 257–269. [11] B. Kalantari, I. Kalantari y R. Zaare-Nahandi, A basic family of iteration functions for polynomial root finding and its characterizations, J. Comput. Appl. Math. 80 (1997), 209–226. [12] D. Kincaid y W. Cheney, An´ alisis num´ erico. Las matem´ aticas del c´ alculo cient´ıfico, AddisonWesley Iberoamericana, Wilmington, 1994. [13] A. Melman, Geometry and convergence of Euler’s and Halley’s methods, SIAM Rev. 39 (1997), 728–735.

220

´ ´ J. A. EZQUERRO, J. M. GUTIERREZ, M. A. HERNANDEZ Y M. A. SALANOVA

[14] B. Neta, Several new methods for solving equations, Intern. J. Computer Math. 22 (1988), 265–282. [15] A. M. Ostrowski, Solution of equations in euclidean and Banach spaces, Academic Press, Nueva York, 1973. [16] G. S. Salehov, On the convergence of the process of tangent hyperbolas (en ruso), Doklady Akad. Nauk SSSR (N.S.) 82 (1952), 525–528. [17] T. R. Scavo y J. B. Thoo, On the geometry of Halley’s method, Amer. Math. Monthly 102 (1995), 417–426. [18] J. K. Stewart, Another variation of Newton’s method, Amer. Math. Monthly 58 (1951), 331– 334. [19] J. F. Traub, Iterative methods for solution of equations, Prentice-Hall, Nueva Jersey, 1964. [20] T. J. Ypma, Historical development of the Newton-Raphson method, SIAM Rev. 37 (1995), 531–551. ´ticas y Computacio ´ n, Universidad de La Rioja, Edificio J. L. Vives, Dpto. de Matema ˜o, Spain 26004 Logron Correo electr´ onico: [email protected], [email protected], [email protected], [email protected]