En muchos problemas de ingeniería, inteligencia artificial u otras disciplinas afines tenemos que resolver ecuaciones del tipo \(f(x)=0\), donde \(f\) es una función que dada una cantidad \(x\), nos devuelve \(f(x)\).

Al valor \(\hat{x}\) que cumple que \(f(\hat{x})=0\) se le llama solución de la ecuación, raíz de la función o cero de la misma.

En dichos problemas, saber explícitamente la función \(f\) en muchos casos no es posible, sólo tenemos un algoritmo que dado un valor \(x\), nos devuelve \(f(x)\).

Entonces, dependiendo de nuestro conocimiento de la función \(f\), podremos aplicar un método numérico u otro.

Todos los métodos numéricos que hallan aproximaciones de ceros construyen una sucesión \((x_n)_n\) que queremos que converja hacia el cero \(\hat{x}\) de la función \(f\).

Cuanto mayor sea la velocidad de convergencia de la sucesión \((x_n)_n\), mejor será el método usado.

Vamos a empezar por el “peor” método de todos por lo que hace a la velocidad de convergencia. Sin embargo, dicha convergencia está asegurada.

El método está basado en el Teorema de Bolzano donde recordemos que dice que si la función \(f\) es continua y tenemos dos valores \(a\) y \(b\) (\(a<b\)) tal que \(f(a)\cdot f(b)<0\), o sea, hay un cambio de signo entre \(a\) y \(b\) o en el intervalo \((a,b)\), podemos asegurar que existe un valor \(\hat{x}\) tal que \(f(\hat{x})=0\).

El método, también llamado método de los intervalos encajados, va construyendo una sucesión de intervalos encajados: \[ [a_0,b_0] \supset [a_1,b_1] \supset \cdots \supset [a_n,b_n]\supset \cdots, \] tal que el cero \(\hat{x}\) siempre está en todos los intervalos \([a_n,b_n]\) y la longitud de cada intervalo \([a_n,b_n]\) vale \(\frac{b_0-a_0}{2^n}\). De esta manera el cero \(\hat{x}\) se calcula de forma más precisa.

¿Cómo se construyen los intervalos \([a_n,b_n]\)?

El primer intervalo \([a_0,b_0]\) vale \([a_0,b_0]=[a,b]\).

A continuación, sea \(c=\frac{a_0+b_0}{2}\) el punto medio del intervalo \([a_0,b_0]\).

Si \(f(a_0)\cdot f(c)<0\), consideramos \(b_1=c\) y \([a_1,b_1]=[a_0,c]\) y si \(f(b_0)\cdot f(c)<0\), consideramos \(a_1=c\) y \([a_1,b_1]=[c,b_0]\).

En general, veamos cómo construir \([a_{n+1},b_{n+1}]\) en función del intervalo \([a_n,b_n]\). Fijaos que se cumplirá siempre que \(f(a_n)\cdot f(b_n)<0\).

Sea \(c=\frac{a_n+b_n}{2}\) el punto medio del intervalo \([a_n,b_n]\). Si \(f(a_n)\cdot f(c)<0\), consideramos \(b_{n+1}=c\) y \([a_{n+1},b_{n+1}]=[a_n,c]\) y si \(f(b_n)\cdot f(c)<0\), consideramos \(a_{n+1}=c\) y \([a_{n+1},b_{n+1}]=[c,b_n]\).

• INPUT a, b, TOL. (damos los valores \(a\) y \(b\) con \(f(a)\cdot f(b)<0\) y la tolerancia admitida donde el cero aproximado \(\hat{x}\) tiene que cumplir \(|f(\hat{x})|\leq \mathrm{TOL}\))
  
• While |f((a+b)/2)|>=TOL do (mientras el punto medio del intervalo en cuestión no cumpla la condición de cero aproximado…)
  • Set c=(a+b)/2. (consideramos el punto medio del intervalo)
  • If f(a)*f(c)<0 (si \(f(a)\cdot f(c)<0\), hemos que definir el nuevo \(b\) como \(c\))
    • Set b=c.
      
  • Else (en caso contrario, el nuevo \(a\) será \(c\))
    • Set a=c.
• Set x=c. (damos como valor aproximado del cero el último valor de \(c\) hallado)

La sucesión \((x_n)_n\) del método de la bisección serán los puntos medios de los intervalos encajados, \(x_n=\frac{a_n+b_n}{2}\).

Como \(x_n,\hat{x}\in [a_n,b_n]\) tendremos que \(|x_n-\hat{x}|\leq b_n-a_n =\frac{b-a}{2^n}\).

Por tanto el orden de convergencia de la sucesión \((x_n)_n\) será la de la sucesión \(\frac{1}{2^n}\): \[ |x_n-\hat{x}|\leq K\cdot \frac{1}{2^n}, \] con \(K=b-a\).

Vamos a considerar otro pseudocódigo para el método cambiando el criterio de parada. Consideraremos que si \(|x_{n+1}-x_n|\leq \mathrm{TOL}\), entonces ya hemos llegado a la aproximación deseada.

• INPUT a, b, TOL. (damos los valores \(a\) y \(b\) con \(f(a)\cdot f(b)<0\) y la tolerancia admitida donde la sucesión que tiende al cero aproximado \(\hat{x}\) tiene que cumplir \(|x_{n+1}-x_n|\leq \mathrm{TOL}\)).
• Set x=(a+b)/2. (definimos el valor \(x_0\))
• Set y=x+1. (el valor \(y\) será donde guardaremos el siguiente elemento de la sucesión. Añadiendo un \(1\), nos aseguramos que “entre” en el bucle While)

• While abs(x-y)>=TOL do (mientras el punto medio del intervalo en cuestión no cumpla la condición de cero aproximado…)
  • Set x=(a+b)/2. (consideramos el punto medio del intervalo)
  • If f(a)*f(x)<0 (si \(f(a)\cdot f(x)<0\), hemos que definir el nuevo \(b\) como \(x\))
    • Set b=x.
  • Else (en caso contrario, el nuevo \(a\) será \(x\))
    • Set a=x.
  • Set y=(a+b)/2. (definimos el nuevo miembro de la sucesión)
• Set cero=x. (damos como valor aproximado del cero el último valor de \(x\) hallado)

¿Cuál de los dos métodos de parada es el mejor?

Los dos métodos son equivalentes pero en el caso en que evaluar la función \(f\) sea muy costoso es mejor usar el segundo método de parada ya que nos evita una evaluación de la función \(f\).

En cambio, si evaluar la función \(f\) no es costoso, podemos usar cualquiera de los dos métodos indistintamente.

El método del punto fijo consiste en transformar la ecuación \(f(x)=0\) en la ecuación \(x=g(x)\) mediante operaciones algebraicas “básicas”.

Entonces el cero de la ecuación \(f(x)=0\) se transforma en lo que llamamos un punto fijo de la función \(g(x)\):

Definición de punto fijo. Sea \(g\) una función real de variable real. Diremos que \(\hat{x}\) es un punto fijo de la función \(g\) si \(g(\hat{x})=\hat{x}\).

Observación: Gráficamente, un punto fijo resulta de la intersección de la recta diagonal \(\color{red}{y=x}\) y de la función \(\color{red}{y=g(x)}\).

La idea es considerar la sucesión definida de forma recurrente como \(x_n = g(x_{n-1})\) y ver bajo qué condiciones la sucesión \((x_n)_n\) converge hacia el punto fijo \(\hat{x}\) de \(g\) que recordemos será el cero de la función \(f\) buscado.

El teorema siguiente nos dice cuáles son las condiciones de existencia del punto fijo de una función \(g\):

Teorema.

• Sea \(g\in {\cal C}([a,b])\) una función continua dentro de un intervalo \([a,b]\) tal que \(g(x)\in [a,b]\), es decir, la función \(g\) está definida de la forma siguiente: \[ g: [a,b]\longrightarrow [a,b]. \] Entonces \(g\) tiene al menos un punto fijo \(\hat{x}\) en el intervalo \([a,b]\), es decir, existe \(\hat{x}\in [a,b]\) tal que \(g(\hat{x})=\hat{x}.\)

• Si además \(g\in {\cal C}^1([a,b])\) tal que existe una constante \(k<1\) con \(|g'(x)|\leq k\), para todo valor \(x\in (a,b)\). Entonces el punto fijo en el intervalo \([a,b]\) es único, es decir, existe un único \(\hat{x}\in [a,b]\) tal que \(g(\hat{x})=\hat{x}.\)

Para resolver la ecuación \(x=g(x)\) o para hallar un punto fijo de la función \(g\), definimos la sucesión \(x_n =g(x_{n-1})\). El Teorema siguiente nos dice bajo qué condiciones la sucesión anterior converge hacia el punto fijo \(\hat{x}\): \(\displaystyle \lim_{n\to\infty}x_n =\hat{x}\).

Teorema del punto fijo. Sea \(g\in {\cal C}([a,b])\) una función continua dentro de un intervalo \([a,b]\) tal que \(g(x)\in [a,b]\). Supongamos, ademas, que existe \(g'(x)\), para todo \(x\in [a,b]\) y que existe una constante \(k\), \(0<k<1\), tal que \(\displaystyle |g'(x)|\leq k\), para todo \(x\in [a,b]\).

Sea \(x_0\in [a,b]\). Definimos la sucesión \((x_n)_n\) de forma recurrente como \(x_n=g(x_{n-1})\), \(n\geq 1\). Entonces \(x_n\) converge hacia el único punto fijo \(\hat{x}\) de \(g(x)\). Es decir \(\displaystyle \lim_{n\to\infty}x_n =\hat{x}\), con \(g(\hat{x})=\hat{x}\).

• INPUT x0, TOL, N. (damos el valor inicial, la tolerancia permitida y el número máximo de iteraciones)
  
• Set iter=1. (iter será un contador que nos dirá el valor de \(n\) en la sucesión \(x_n\))

• While iter <= N do (mientras no llegemos al máximo de iteraciones hacer lo siguiente)
  • Set x=g(x0). (calculamos el siguiente término de la sucesión \(x_n\))
    
  • If |x-x0| < TOL then (miramos si el valor \(x_n\) cumple la condición de la tolerancia permitida)
    • print (x); (damos la solución)
      
    • stop. (acabamos)
      
  • iter++. (aumentamos el valor \(n\) de la sucesión)
    
  • Set x0=x (actualizamos el valor \(x_n\))
    
• Print('El método ha fallado después de N iteraciones') (si hemos llegado hasta aquí significa que hemos alcanzado el número máximo de iteraciones sin que se cumpla la condición de la tolerancia)

En la demostración del Teorema del punto fijo vimos una estimación del error cometido con la aproximación \(x_n\) con respecto al punto fijo \(\hat{x}\): \[ |x_n-\hat{x}|\leq k^n\cdot |x_0-\hat{x}|. \] Dicha estimación no es útil ya que al no conocer \(\hat{x}\), tampoco conocemos \(|x_0-\hat{x}|.\)

La proposición siguiente nos da una estimación del error \(|x_n-\hat{x}|\) en función de los extremos del intervalo \([a,b]\), \(x_0\) y de \(x_1\) que son valores conocidos.

Proposición. Estimación del error cometido. En las condiciones del Teorema del punto fijo, se cumple: \[ |x_n-\hat{x}|\leq k^n\cdot\max\{x_0-a,b-x_0\},\quad |x_n-\hat{x}|\leq \frac{k^n}{1-k}\cdot |x_1-x_0|, \] para todo \(n\geq 1\), si \(x_0\in [a,b]\).

El método de Newton-Raphson es uno de los métodos más conocidos y más potentes en el sentido de su velocidad de convergencia. Sin embargo, al contrario que el método de la bisección, su convergencia no está asegurada.

Veamos gráficamente en qué consiste el método de Newton-Raphson:

En el gráfico anterior vemos que tratamos de hallar un cero \(\hat{x}\) marcado en rojo de la función \(f(x)=\mathrm{e}^{-x}-\frac{1}{2}\ln x\).

Para ello, empezamos con el valor inicial \(x_0=1\). A continuación hallamos la recta tangente a la función \(f(x)\) en el punto \((x_0,f(x_0))=\left(1,\mathrm{e}^{-1}\right)\approx (1,0.367879)\).

Seguidamente, miramos dónde corta dicha recta tangente al eje de las \(X\). Dicho punto de corte, será el siguiente elemento de la sucesión \(x_1\).

El paso siguiente es hacer con el punto \(x_1\) lo mismo que hemos hecho con el punto \(x_0\): hallar la recta tangente a la función \(f(x)\) en el punto \((x_1,f(x_1))\) y mirar dónde corta dicha recta tangente al eje de las \(X\) para hallar \(x_2\) y así sucesivamente.

De esta forma obtenemos una sucesión \(x_0,x_1,x_2,\ldots\) cuyo límite como puede observarse en el gráfico es el cero \(\hat{x}\).

Como puede verse en el gráfico, la velocidad de convergencia es muy rápida ya que con dos iteraciones, tenemos un valor \(x_2\) muy próximo al cero \(\hat{x}\).

El método de Newton-Raphson que nos da el valor de \(x_{n}\) en función de \(x_{n-1}\) se deduce de realizar los pasos siguientes:

• Desarrollamos por Taylor la función \(f\) que suponemos de clase \({\cal C}^2[a,b]\) en un cierto intervalo de la que queremos hallar el cero en el punto \(x_{n-1}\) usando el polinomio de Taylor de grado \(1\) y usando la expresión del resto de Lagrange: \[ f(x)=f(x_{n-1})+f'(x_{n-1})\cdot (x-x_{n-1})+\frac{f''(\xi(x_{n-1}))}{2}\cdot (x-x_{n-1})^2, \] donde \(\xi(x_{n-1})\in <x,x_{n-1}>\).

• Si usamos la aproximación anterior en el cero \(\hat{x}\), como \(f(\hat{x})=0\), tenemos que: \[ 0=f(\hat{x})=f(x_{n-1})+f'(x_{n-1})\cdot (\hat{x}-x_{n-1})+\frac{f''(\xi(x_{n-1}))}{2}\cdot (\hat{x}-x_{n-1})^2. \]

• Si suponemos que \(|\hat{x}-x_{n-1}|\) es pequeño podemos despreciar el término del resto de Lagrange y obtenemos: \[ 0=f(x_{n-1})+f'(x_{n-1})\cdot (x-x_{n-1}), \] donde hemos escrito \(x\) en lugar de \(\hat{x}\) ya que el valor que cumple la condición anterior no és el cero \(\hat{x}\) sino un valor aproximado \(x\). Dicho valor aproximado será el siguiente valor de la sucesión \(x_n\):

• Si despejamos \(x=x_n\) de la ecuación anterior, hallamos la relación buscada: \[ x_n =x_{n-1}-\frac{f(x_{n-1})}{f'(x_{n-1})}. \]

Observación: Los pasos anteriores equivalen a buscar la recta tangente en el punto \((x_{n-1},f(x_{n-1}))\), cuya ecuación es \(y=f(x_{n-1})+f'(x_{n-1})\cdot (x-x_{n-1})\), y buscar el corte de dicha recta con el eje \(X\), es decir, resolver \(y=0\), para hallar el punto \(x_n\), tal como comentamos en la introducción del método.

Observación: El método de Newton-Raphson es un caso particular de usar el metódo del punto fijo, con la función \(g(x)=x-\frac{f(x)}{f'(x)}\). Observar que el valor \(x_n\) se calcular de la forma siguiente en función de la función \(g\): \(x_n=g(x_{n-1})\).

• INPUT x0, TOL, N. (damos el valor inicial \(x_0\), la tolerancia permitida y el número máximo de iteraciones \(N\))
  
• Set n=1. (\(n\) será el contador para las iteraciones realizadas)

• While n<=N do (mientras no llegemos al número máximo de interaciones hacer lo siguiente)
  • Set x=x0-f(x0)/f'(x0). (calculamos el siguiente valor de la iteración)
    
  • If |x-x0|< TOL then (si la diferencia entre los dos valores verifica la condición de la tolerancia…)
    • Print(x); (dar la aproximación del cero)
      
    • STOP. (acabar)
      
  • Set n++. (si no acabamos aumentar el contador del número de iteraciones en \(1\))
    
  • Set x0=x (actualizamos el valor de x0)
• Print ('el método ha fallado después de N iteraciones) (si hemos llegado hasta aquí significa que el método no converge en nuestro caso)
  
• STOP.

Ya hemos comentado que no tenemos asegurada la convergencia del método de Newton-Raphson. Sin embargo, el Teorema siguiente nos dice bajo qué condiciones podemos tener convergencia:

Teorema. Sea \(f\in {\cal C}^2([a,b])\). Sea \(\hat{x}\) un cero de la función \(f\), \(f(\hat{x})=0\) con \(f'(\hat{x})\neq 0\), es decir \(\hat{x}\) es un cero simple de \(f\). Entonces existe un valor \(\delta >0\) tal que la sucesión generada por el método de Newton-Raphson, \((x_n)_n\) converge a \(\hat{x}\) para todo valor inicial \(x_0\in (\hat{x}-\delta,\hat{x}+\delta)\).

El teorema anterior nos dice que el método de Newton-Raphson converge siempre que el cero \(\hat{x}\) sea simple, es decir, \(f'(\hat{x})\neq 0\).

La velocidad de convergencia dependerá de la cota \(k\) de la función \(g'(x)=\frac{f(x)\cdot f''(x)}{f'(x)^2}\) en un entorno de radio \(\delta\) del cero \(\hat{x}\). Dicho valor \(k\) en general será desconocido y hallarlo es un problema mucho más difícil que hallar el cero \(\hat{x}\).

Por tanto, en la práctica, se realiza un gráfico aproximado de la función \(f(x)\) y se elige un punto \(x_0\) cerca del valor de corte de la función \(f\) con el eje \(X\).

El método de Newton es un método muy potente pero tiene un inconveniente: necesitamos conocer la función derivada \(f'(x)\).

En muchos problemas de análisis numérico, dicho conocimiento es un lujo. Sabemos cómo evaluar la función \(f\) pero no hay forma de tener una expresión de la función derivada \(f'\). En estos casos, el metodo de Newton-Raphson no es aplicable.

Para solventar esta dificultad, podemos usar el cociente incremental como aproximación a \(f'(x_{n-1})\): \[ f'(x_{n-1})\approx \frac{f(x_{n-1})-f(x_{n-2})}{x_{n-1}-x_{n-2}}, \] si \(x_{n-2}\) y \(x_{n-1}\) están próximos, la aproximación anterior tendrá poco error.

Si sustituimos la expresión anterior de \(f'(x_{n-1})\) en el método de Newton-Raphson, obtenemos el denominado metodo de la secante: \[ x_n =x_{n-1}-\frac{f(x_n)}{\frac{f(x_{n-1})-f(x_{n-2})}{x_{n-1}-x_{n-2}}}=x_{n-1}-\frac{f(x_{n-1})\cdot (x_{n-1}-x_{n-2})}{f(x_{n-1})-f(x_{n-2})}. \]

El método de la secante también puede deducirse siguiendo los pasos siguientes:

• Dados los puntos \((x_{n-2},f(x_{n-2}))\) y \((x_{n-1},f(x_{n-1}))\), consideramos la recta que pasa por dichos puntos de ecuación: \[ y-f(x_{n-1})=\frac{f(x_{n-1})-f(x_{n-2})}{x_{n-1}-x_{n-2}}\cdot (x-x_{n-1}). \]

• Hallamos la intersección de dicha recta con el eje \(X\) resolviendo \(y=0\) y el punto de intersección será el punto siguiente de la sucesión \(x_n\): \[ y=0,\ \Rightarrow x=x_n=x_{n-1}-\frac{f(x_{n-1})\cdot (x_{n-1}-x_{n-2})}{f(x_{n-1})-f(x_{n-2})}, \] obteniendo la misma expresión que anteriormente.

En el gráfico anterior vemos cómo funciona el método de la secante:

• Hemos considerado \(x_0=0.5\) y \(x_1=1\). Hallando la recta que pasa por los puntos \((0.5,f(0.5))\) y \((1,f(1))\) (en rojo de puntitos) y viendo dónde corta el eje de las \(X\), hallamos el nuevo punto de la sucesión \(x_2\).

• Hacemos lo mismo con los puntos \((x_1,f(x_1))=(1,f(1))\) y \((x_2,f(x_2))\), es decir, hallamos la recta que pasa por los puntos \((1,f(1))\) y \((x_2,f(x_2))\) (en rojo de puntitos) y viendo dónde corta el eje de las \(X\), hallamos el nuevo punto de la sucesión \(x_3\).

• Y asi sucesivamente, hasta llegar a una aproximación del cero \(\hat{x}\).

• INPUT x0, x1, TOL, N. (Damos los valores iniciales \(x_0\), \(x_1\), la tolerancia TOL y el número máximo de iteraciones \(N\))
  
• Set n=2. (El valor \(n\) nos da el subíndice del nuevo \(x_n\) a calcular)
  
• Set y0=f(x0). (Definimos \(y_0=f(x_0)\))
  
• Set y1=f(x1). (Definimos \(y_0=f(x_1)\))

• While n <= N do (Mientras no lleguemos al número de máximo de iteraciones hacer lo siguiente…)
  • Set x=x1-y1*(x1-x0)/(y1-y0). (Calculamos el nuevo valor \(x_n\))
    
  • If |x-x1| < TOL then (Si el valor \(x_n\) cumple la condición de la tolerancia…)
    • Print(x) (Damos \(x_n\) como aproximación del cero \(\hat{x}\))
      
    • STOP (Acabamos)
      
  • Set n++. (En caso contrario, aumentamos \(n\) en \(1\))
    
  • Set x0=x1. (Actualizamos \(x_0, x_1, y_0\) e \(y_1\))
    
  • Set y0=y1.
    
  • Set x1=x.
  • Set y1=f(x).

• Print ('El método ha fallado después de N iteraciones). (Si hemos llegado hasta aquí significa que el método ha fallado y ha alcanzado el número máximo de iteraciones)
  
• STOP. (Paramos)

Los métodos de Newton-Raphson y de la secante son método relativamente rápidos en cuanto a su convergencia pero tienen el inconveniente que no tenemos asegurada dicha convergencia.

Recordemos el Teorema donde nos dice si el cero es simple (\(f'(\hat{x})\neq 0\)), existe un entorno del cero \(\hat{x}\) tal que para cualquier valor inicial \(x_0\) perteneciente a dicho entorno, el método de Newton-Raphson es convergente. El problema es que hallar dicho entorno es un problema mucho más difícil que hallar el cero en sí.

Por dicho motivo, de cara a asegurar la convergencia de la sucesión \((x_n)_n\) tenemos el método de regula falsi que usa las ventajas de los métodos de la bisección y de la secante.

En pocas palabras, el método de regula falsi supone que hay un cambio de signo entre \(x_{n-2}\) y \(x_{n-1}\), es decir, \(f(x_{n-2})\cdot f(x_{n-1}) <0\) y halla \(x_n\) usando el método de la secante. Como \(x_n\) estará entre \(x_{n-2}\) y \(x_{n-1}\) habrá un cambio de signo de \(f\) entre \(x_{n-2}\) y \(x_n\) o entre \(x_{n-1}\) y \(x_n\). Dependiendo de dónde esté el cambio de signo, cambia el valor de \(x_{n-1}\).

Más concretamente, el método de regula falsi consiste en los pasos siguientes: (vamos a dar los pasos para pasar al paso \(n\) desde el paso \(n-1\))

• Sean \(x_{n-2}\) y \(x_{n-1}\) tal que \(f(x_{n-2})\cdot f(x_{n-1})<0\). No suponemos ningún orden sobre \(x_{n-2}\) y \(x_{n-1}\), es decir, podemos tener \(x_{n-2}<x_{n-1}\) o \(x_{n-1}<x_{n-2}\).

• Calculamos \(x_n\) según el método de la secante: \[ x_n=x_{n-1}-\frac{f(x_{n-1})\cdot (x_{n-1}-x_{n-2})}{f(x_{n-1})-f(x_{n-2})}. \]

• Caso 1. Si \(f(x_n)\cdot f(x_{n-1})<0\), tenemos un cambio de signo de \(f\) en el intervalo \(<x_{n-1},x_n>\). En este caso, no hacemos nada y la sucesión quedará \(\ldots,x_{n-2},x_{n-1},x_n,\ldots\)
• Caso 2. Si \(f(x_n)\cdot f(x_{n-2})<0\) tenemos un cambio de signo de \(f\) en el intervalo \(<x_{n-2},x_n>\). En este caso, cambiamos el orden de \(x_{n-2}\) y \(x_{n-1}\) y la sucesión quedará \(\ldots,x_{n-1},x_{n-2},x_n,\ldots\)
• Volvemos al paso inicial para hallar \(x_{n+1}\).
• Y así sucesivamente.

El gráfico anterior ilustra el método de regula falsi.

• Empezamos con \(x_0=1\) y \(x_1=2\) donde observamos que hay un cambio de signo de la función \(f\). Más concretamente, \(f(x_0)<0\) y \(f(x_1)>0\).

• Hallamos el punto \(x_2\) usando el método de la secante y vemos que estamos en el caso 2. La sucesión será por tanto \(x_1=2,x_0=1,x_2,\ldots\)

• Hallamos \(x_3\) usando el método de la secante con los puntos \((x_0,f(x_0))=(1,f(1))\) y \((x_2,f(x_2))\) y vemos que volvemos a estar en el caso 2. La sucesión será por tanto \(x_1=2,x_2,x_0=1,x_3\).

• Y así sucesivamente.

Observamos que de esta forma la sucesión \((x_n)_n\) siempre converge hacia el cero \(\hat{x}\), hecho que no teníamos asegurado con el método de la secante.

• INPUT x0, x1, TOL, N. (Damos los valores iniciales x0 y x1, la tolerancia TOL y el número máximo de iteraciones \(N\))
  
• Set n=2. (El valor \(n\) nos da el subíndice del nuevo \(x_n\) a calcular)
  
• Set y0=f(x0). (Definimos \(y_0=f(x_0)\))
  
• Set y1=f(x1). (Definimos \(y_1=f(x_1)\))

• While n <= N do (Mientras no lleguemos al número de máximo de iteraciones hacer lo siguiente…)
  • Set x=x1-y1*(x1-x0)/(y1-y0). (Calculamos el nuevo valor \(x_2\))
    • If min(|x-x1|, |x-x0|) < TOL then (Si el valor \(x_2\) cumple la condición de la tolerancia…)
      
    • Print(x) (Damos \(x_2\) como aproximación del cero \(\hat{x}\))
      
    • STOP (Acabamos)
      
  • Set n++. (En caso contrario, aumentamos \(n\) en \(1\))
  • Set y=f(x). (Definimos \(y_2=f(x_2)\))

• • If y*y1<0 then (Si el cambio de signo está en el intervalo \(<x_1,x_2>\))
    • Set x0=x1. (El nuevo \(x_0\) será el “viejo” \(x_1\))
      
    • Set y0=y1.
  • Set x1=x. (El nuevo \(x_1\) será el valor de \(x_2\))
    
  • Set y1=y.
    
• Print ('El método ha fallado después de N iteraciones). (Si hemos llegado hasta aquí significa que el método ha fallado y ha alcanzado el número máximo de iteraciones)
  
• STOP. (Paramos)

El pseudocódigo necesita un poco de explicación respecto del algoritmo que hemos dados antes.

Básicamente, empezamos con dos valores x0 y x1 donde hay un cambio de signo de la \(f\) y calculamos el nuevo valor x.

A continuación, tenemos que redefinir los nuevos x0 y x1 dependiende de dónde esté el cambio de signo:

• Si estamos en el caso 1, es decir, \(f(\mathtt{x1})\cdot f(\mathtt{x})<0\), basta que definamos el “nuevo” x0 como x1 y el “nuevo” x1 como x y tendremos que volverá a haber un cambio de signo entre x0 y x1.
• Si estamos en el caso 2, es decir, \(f(\mathtt{x0})\cdot f(\mathtt{x})<0\), basta que definamos el “nuevo” x1 como x y tendremos que volverá a haber un cambio de signo entre x0 y x1.

Observamos que la convergencia es más lenta que el método de la secante debido básicamente a que tenemos que asegurarnos del cambio de signo en el intervalo \(<x_{n-2},x_{n-1}>\).

Sin embargo, la convergencia está asegurada siempre, cosa que no podemos decir lo mismo del método de la secante.

Cuando un método numérico consiste en hallar una sucesión \((x_n)_n\) que converge hacia una solución \(\hat{x}\), un concepto importante para medir lo rápido de converge es el orden de convergencia de la sucesión:

Definición de orden de convergencia. Sea \((x_n)_n\) una sucesión de números reales convergente hacia \(\hat{x}\), donde suponemos que \(x_n\neq \hat{x}\), para cualquier \(n\). Suponemos que existen constantes \(\lambda >0\) y \(\alpha>0\) tal que \[ \lim_{n\to\infty} \frac{|x_{n+1}-\hat{x}|}{|x_n-\hat{x}|^\alpha} =\lambda. \] Entonces, diremos que \((x_n)_n\) converge hacia \(\hat{x}\) con orden \(\alpha\) con constante asintótica \(\lambda\).

• Cuando mayor es el orden \(\alpha\), más rápido converge la sucesión \((x_n)_n\) hacia \(\hat{x}\). La idea es que para \(n\) grande, podemos aproximar el error en el término \(n+1\) en función del error en el término \(n\) de la forma siguiente: \[ |x_{n+1}-\hat{x}|\approx \lambda |x_{n}-\hat{x}|^\alpha. \] Entonces para valores grandes de \(\alpha\) el error irá disminuyendo muy rápidamente, acercándose la sucesión \((x_n)_n\) a su límite \(\hat{x}\) de forma muy rápida.

• Si \(\alpha=1\) y \(\lambda\leq 1\), se dice que la sucesión \((x_n)_n\) converge linealmente.

• Si \(\alpha=2\), se dice que la sucesión \((x_n)_n\) converge de forma cuadrática.

• El orden de convergencia de la mayoría de las sucesiones que “tenemos en mente” tienen convergencia lineal:
  • \(x_n = \frac{P(n)}{Q(n)}\), donde \(P\) y \(Q\) son dos polinomios del mismo grado en \(n\), donde \(\hat{x}=\frac{a}{b}\) siendo \(a\) y \(b\) los coeficientes de mayor grado de \(P\) y \(Q\), respectivamente y \(\lambda=1\),
  • \(x_n= a^n\), con \(a<1\) donde \(\hat{x}=0\) y \(\lambda =a\),
  • \(x_n= \left(1+\frac{1}{n}\right)^n\), donde \(\hat{x}=\mathrm{e}\) y \(\lambda=1\).
• Los casos interesantes son aquellos en los que la sucesión de define por recurrencia como es el caso del método del punto fijo: \(x_n=g(x_{n-1})\).

Teorema. Sea \(g\in {\cal C}[a,b]\) una función continua en el intervalo \([a,b]\) tal que \(g(x)\in [a,b]\) para todo valor \(x\in [a,b]\).

Supongamos además que \(g\in {\cal C}^1 (a,b)\) tiene derivada continua en el intervalo \((a,b)\) y que existe una constante \(k<1\) tal que \(|g'(x)|\leq k\), para todo valor \(x\in (a,b)\). Sea \(x_0\in [a,b]\), definimos la sucesión de forma recurrente \(x_n=g(x_{n-1})\) que sabemos que es convergente hacia \(\hat{x}\in [a,b]\).

Entonces, si \(g'(\hat{x})\neq 0\), la sucesión anterior \((x_n)_n\) converge de forma lineal hacia \(\hat{x}\).

El teorema siguiente nos da las condiciones para que el método del punto fijo tenga convergencia cuadrática:

Teorema. Sea \(g\in {\cal C}[a,b]\) una función continua en el intervalo \([a,b]\).

Sea \(\hat{x}\) una solución de la ecuación \(x=g(x)\) tal que \(g'(\hat{x})=0\). Supongamos además que \(g\in {\cal C}^2 (a,b)\) tiene las dos primeras derivadas continuas en el intervalo \((a,b)\) y que existe una constante \(M\) tal que \(|g''(x)|< M\) en un intervalo abierto \(I\subseteq (a,b)\) tal que \(\hat{x}\in I\).

Entonces existe un valor \(\delta >0\) tal que si \(x_0\in (\hat{x}-\delta,\hat{x}+\delta)\), la sucesión definida de forma recurrente por \(x_n=g(x_{n-1})\) converge en orden cuadrático hacia \(\hat{x}\). Además, para \(n\) suficientemente grande: \[ |x_{n+1}-\hat{x}|< \frac{M}{2}\cdot |x_n-\hat{x}|^2. \]

Veamos un resumen de la sección.

Sea \(g\in {\cal C}[a,b]\) una función continua en el intervalo \([a,b]\) tal que tiene las dos primeras derivadas continuas en \((a,b)\). Sea \(x_0\in [a,b]\) y definimos \(x_n=g(x_{n-1})\). Supongamos que \((x_n)_n\) es convergente hacia \(\hat{x}\in [a,b]\). Entonces:

• Si \(g'(\hat{x})\neq 0\) con \(|g'(\hat{x})|<1\), la sucesión \((x_n)_n\) converge linealmente hacia \(\hat{x}\).
• Si \(g'(\hat{x})=0\) y \(g''(\hat{x})\neq 0\), la sucesión \((x_n)_n\) converge cuadráticamente hacia \(\hat{x}\).

Veamos en qué condiciones el método de Newton-Raphson converge cuadráticamente hacia el cero:

Teorema. Sea \(f\in {\cal C}^2([a,b])\). Sea \(\hat{x}\) un cero de la función \(f\), \(f(\hat{x})=0\) con \(f'(\hat{x})\neq 0\), es decir \(\hat{x}\) es un cero simple de \(f\). Entonces la sucesión \((x_n)_n\) generada por el método de Newton-Raphson converge al menos cuadráticamente hacia \(\hat{x}\).

Si nos dan un método iterativo para hallar ceros, como por ejemplo el de la secante o el de la regula-falsi, ¿existe una manera aproximada de calcular el orden de convergencia \(\alpha\) del método anterior?

Supongamos que \(x_n=g(x_{n-1})\). Si el orden de convergencia es \(\alpha\), tendremos que: \[ |x_{n+1}-\hat{x}|\approx \lambda\cdot |x_n-\hat{x}|^\alpha. \] Tomando logaritmos en la expresión anterior obtenemos: \[ \ln |x_{n+1}-\hat{x}|\approx \ln \lambda +\alpha \ln |x_{n}-\hat{x}|. \]

Vemos por tanto que la sucesión \(\ln |x_{n+1}-\hat{x}|\) depende linealmente de la sucesión \(\ln |x_{n}-\hat{x}|\).

Es decir, que si hallamos la recta de regresión de la sucesión \(\ln |x_{n+1}-\hat{x}|\) como función de la sucesión \(\ln |x_{n}-\hat{x}|\) la pendiente de dicha recta será aproximadamente el valor del orden de convergencia \(\alpha\).

El método anterior tiene una pega: ¡no conocemos el cero \(\hat{x}\)! Ahora bien, como \(\displaystyle\lim_{n\to\infty}x_n =\hat{x}\), podemos aproximar \(\hat{x}\approx x_N\) con \(N\) valor grande.

Una raíz múltiple de una función \(f\) es aquélla en la que la función y las derivadas hasta cierto orden son cero:

Definición de raíz múltiple. Sea \(f\in {\cal C}^{k}[a,b]\) una función de clase \(k\) que significa que la función y las \(k\) primeras derivadas son continuas. Diremos que \(\hat{x}\in [a,b]\) es un cero múltiple de multiplicidad \(k\) si \(f(\hat{x})=f'(\hat{x})=\ldots =f^{(k-1)}(\hat{x})=0\) y \(f^{(k)}(\hat{x})\neq 0\).

Los ceros de multiplicidad \(1\) se denominan ceros simples, son los que verifican que \(f(\hat{x})=0\) pero \(f'(\hat{x})\neq 0\).

Los ceros de multiplicidad \(2\) se denominan ceros dobles, son los que verifican que \(f(\hat{x})=f'(\hat{x})=0\) pero \(f''(\hat{x})\neq 0\).

Observación. Otra forma de definir la multiplicidad de un cero \(\hat{x}\) es decir que \(\hat{x}\) es de multiplicidad k si, y sólo si, \(f(x)=(x-\hat{x})^k q(x)\), si \(x\neq \hat{x}\), con \(q(\hat{x})\neq 0\).

La demostración de la equivalencia de las definiciones se basa en desarrollar por Taylor la función \(f\) en un entorno de \(\hat{x}\).

La velocidad de convergencia de un método depende de la multiplicidad del cero en cuestión. Por ejemplo, para que el método de Newton-Raphson tenga orden de convergencia \(2\) necesitamos que \(\hat{x}\) sea un cero simple tal como nos dice el Teorema correspondiente.

Si observamos que un método de convergencia rápida como por ejemplo el método de Newton-Raphson tiene velocidad lenta a la hora de calcular un cero de una función \(f\), podemos sospechar que el cero es una raíz múltiple de \(f\). En estos casos, podemos considerar la función \(\tilde{f}(x)=\frac{f(x)}{f'(x)}\) en lugar de la función \(f\):

Proposición. Sea \(f\in {\cal C}^k (a,b)\) una función de clase \({\cal C}^k\) en el intervalo \((a,b)\) tal que \(\hat{x}\) es una raíz de multiplicidad \(k\). Entonces \(\hat{x}\) es una raíz simple de la función \(\tilde{f}(x)=\frac{f(x)}{f'(x)}\).

Usando la proposición anterior, el método de Newton-Raphson para calcular ceros múltiples se reduce a hallar una sucesión \((x_n)_n\) que verifica: \[ \begin{align*} x_n= & x_{n-1}-\frac{\tilde{f}(x_{n-1})}{\tilde{f}'(x_{n-1})}=x_{n-1}-\frac{\frac{f(x_{n-1})}{f'(x_{n-1})}}{\frac{f'(x_{n-1})^2-f(x_{n-1})\cdot f''(x_{n-1})}{f'(x_{n-1})^2}}\\ = & x_{n-1}-\frac{f(x_{n-1})\cdot f'(x_{n-1})}{f'(x_{n-1})^2-f(x_{n-1})\cdot f''(x_{n-1})}. \end{align*} \]

Tener convergencia cuadrática en un método para hallar ceros de funciones no siempre es posible ya que puede pasar que el cero sea múltiple o que no sea posible evaluar las derivadas de la función.

Por ejemplo, hay veces que aunque podamos evaluar la derivada primera de la función, no es posible o muy complicado evaluar las derivadas de orden superior.

Recordemos que si queremos aplicar el método de Newton-Raphson para el cálculo de ceros múltiples, es necesario poder evaluar la derivada primera y la segunda, cosas que en muchos casos es un lujo.

Por dichos motivos, vamos a ver dos métodos para acelerar la convergencia de sucesiones convergentes: el método de Aitken y el método de Steffensen.

Sea \((x_n)_n\) una sucesión convergente hacia \(\hat{x}\) con convergencia lineal.

Vamos a construir otra sucesión \(y_n\) que converge más rápidamente que \(x_n\) hacia \(\hat{x}\).

Para ello, supongamos que para \(n\) suficientemente grande, \[ \frac{x_{n+1}-\hat{x}}{x_n -\hat{x}}\approx \frac{x_{n+2}-\hat{x}}{x_{n+1} -\hat{x}}. \]

La suposición anterior no es descabellada ya que sabemos que

• como \((x_n)_n\) es convergente, los signos de \(x_{n+1}-\hat{x}\), \(x_{n}-\hat{x}\), \(x_{n+2}-\hat{x}\) coincidirán y
• como tiene convergencia lineal, existe un valor \(\lambda >0\) tal que \[ \lim_{n\to \infty}\frac{|x_{n+1}-\hat{x}|}{|x_n -\hat{x}|}=\lim_{n\to \infty} \frac{x_{n+1}-\hat{x}}{x_n -\hat{x}} = \lim_{n\to \infty} \frac{x_{n+2}-\hat{x}}{x_{n+1} -\hat{x}}= \lambda >0 \]

Por tanto, \[ \begin{align*} (x_{n+1}-\hat{x})^2\approx & (x_{n+2}-\hat{x})\cdot (x_n-\hat{x}),\\ x_{n+1}^2-2x_{n+1}\hat{x}+\hat{x}^2 \approx & x_{n+2}x_n-(x_n+x_{n+2})\hat{x}+\hat{x}^2,\\ (x_{n+2}+x_n-2x_{n+1})\hat{x}\approx & x_{n+2}x_n-x_{n+1}^2,\\ \hat{x}\approx & \frac{x_{n+2}x_n-x_{n+1}^2}{x_{n+2}+x_n-2x_{n+1}} =x_n -\frac{(x_{n+1}-x_n)^2}{x_{n+2}-2x_{n+1}+x_n}, \end{align*} \] donde la última igualdad se demuestra haciendo manipulaciones algebraicas básicas.

Definición del método de Aitken. El método de Aitken consiste en transformar la sucesión convergente \((x_n)_n\) hacia \(\hat{x}\) en otra sucesión \((y_n)_n\) de la forma siguiente: \[ y_n=x_n -\frac{(x_{n+1}-x_n)^2}{x_{n+2}-2x_{n+1}+x_n}, \] cuya convergencia es más rápida que la sucesión original \((x_n)_n\).

Antes de pasar a demostrar que la sucesión \((y_n)_n\) converge más rápidamente que la sucesión \((x_n)_n\) vamos a introducir alguna notación:

Definición de diferencia hacia delante \(\Delta x_n\). Sea \((x_n)_n\) una sucesión de números reales. La sucesión de diferencias hacia delante \((\Delta x_n)_n\) se define como: \[ \Delta x_n =x_{n+1}-x_n, \mbox{ para }n\geq 0, \] y la sucesión de diferencias hacia delante de orden \(k\), para \(k\geq 2\), \((\Delta^k x_n)_n\), se define de forma recurrente como: \[ \Delta^k x_n = \Delta(\Delta^{k-1} x_n),\mbox{ para }n\geq 0. \]

Por ejemplo, la sucesión de diferencias hacia delante de orden \(2\), \((\Delta^2 x_n)_n\) sería: \[ \begin{align*} \Delta^2 x_n = & \Delta(\Delta x_n)=\Delta x_{n+1}-\Delta x_n =(x_{n+2}-x_{n+1})-(x_{n+1}-x_n)\\ = & x_{n+2}-2 x_{n+1}+x_n,\mbox{ para }n\geq 0. \end{align*} \]

Usando la notación introducida anteriormente, la sucesión \((y_n)_n\) del método de Aitken se puede escribir como: \[ y_n = x_n -\frac{(x_{n+1}-x_n)^2}{x_{n+2}-2x_{n+1}+x_n} =x_n-\frac{(\Delta x_n)^2}{\Delta^2 x_n}. \] La expresión anterior es el motivo por el que el método de Aitken se denomina método \(\Delta^2\) de Aitken.

La proposición siguiente nos dice bajo qué condiciones la sucesión definida por el método \(\Delta^2\) de Aitken converge más rápidamente que la sucesión original:

Proposición. Sea \((x_n)_n\) una sucesión convergente hacia \(\hat{x}\) con convergencia lineal tal que \[ \lim_{n\to\infty} \frac{x_{n+1}-\hat{x}}{x_n-\hat{x}}<1. \] Entonces la sucesión \((y_n)_n\) definida por el método \(\Delta^2\) de Aitken converge más rápidamente hacia \(\hat{x}\) en el sentido de que: \[ \lim_{n\to\infty} \frac{y_{n}-\hat{x}}{x_n-\hat{x}}=0. \]

El método de Steffensen se aplica para acelerar la convergencia de sucesiones generadas usando el método del punto fijo.

Es decir, suponemos que tenemos una sucesión \((x_n)_n\) convergente y una función \(g\) tal que \(x_n=g(x_{n-1})\). Entonces \((x_n)_n\) converge hacia un punto fijo \(\hat{x}\) de \(g(x)\): \(g(\hat{x})=\hat{x}\).

El método genera una sucesión nueva \((y_n)_n\) pero no a partir de la original \((x_n)_n\) sino a partir de la función \(g\). Veamos los pasos que sigue:

• Sean \(x_0,x_1=g(x_0)\) y \(x_2=g(x_1)\) los tres primeros términos de la sucesión \((x_n)_n\) a los que llamaremos \(x_0^{(0)}\), \(x_1^{(0)}\) y \(x_2^{(0)}\), respectivamente.

• El primer término \(x_0^{(1)}\) se calcula aplicando el método de Aitken a los tres términos anteriores: \[ x_0^{(1)}=x_0^{(0)}-\frac{(x_{1}^{(0)}-x_0^{(0)})^2}{x_{2}^{(0)}-2x_{1}^{(0)}+x_0^{(0)}}=x_0^{(0)}-\frac{(\Delta x_0^{(0)})^2}{\Delta^2 x_0^{(0)}}. \]

• Calculamos \(x_1^{(1)}=g(x_0^{(1)})\) e \(x_2^{(1)}=g(x_1^{(1)})\).

• El segundo término \(x_0^{(2)}\) se calcula aplicando el método de Aitken a los tres términos anteriores: \[ x_0^{(2)}=x_0^{(1)}-\frac{(x_{1}^{(1)}-x_0^{(1)})^2}{x_{2}^{(1)}-2x_{1}^{(1)}+x_0^{(1)}}=x_0^{(1)}-\frac{(\Delta x_0^{(1)})^2}{\Delta^2 x_0^{(1)}}. \]

• Y así sucesivamente.

Los términos de la nueva sucesión \((y_n)_n\) serán, pues, \(y_n=x_0^{(n)}\).

Observación. Podría pasar que, usando el algoritmo anterior, \(\Delta^2 x_0^{(n)}=0\) para algún \(n\).

En este caso, no podríamos continuar ya que tendríamos un cero en el denominador a la hora de calcular \(y_{n+1}=x_0^{(n+1)}\).

Si esto sucediera, paramos y consideramos \(x_0^{(n)}\) como la mejor aproximación del cero o del punto fijo de \(g\).

El siguiente resultado nos dice que el método de Steffensen es capaz de transformar una sucesión con convergencia lineal a convergencia cuadrática:

Proposición. Sea \(g\in {\cal C}^1(a,b)\) una función de classe \({\cal C}^1\) en un intervalo \((a,b)\). Supongamos que la ecuación \(x=g(x)\) tiene como solución \(\hat{x}\) con \(g'(\hat{x})\neq 0\).

Si existe un \(\delta>0\) tal que \(g\in {\cal C}^3(\hat{x}-\delta,\hat{x}+\delta)\), el método de Steffensen transforma la sucesión \(x_n=g(x_{n-1})\) con \(x_0\in (\hat{x}-\delta,\hat{x}+\delta)\) de límite \(\hat{x}\) en una sucesión \((y_n)_n\) con convergencia cuadrática.

Observación. La sucesión de \((x_n)_n\) tiene convergencia lineal ya que su límite verifica \(g'(\hat{x})\neq 0\).

• INPUT x0, TOL, N. (Damos como valores iniciales el valor \(x_0\), la tolerancia permitida y el número máximo de iteraciones)
  
• Set n=1. (\(n\) será el contador para las iteraciones realizadas)
  
• While n <= N do (mientras no llegemos al máximo de iteraciones hacer lo siguiente)
  • Set x1=g(x0). (calculamos \(x_1\))
    
  • Set x2=g(x1). (calculamos \(x_2\))
    
  • Set x=x0-(x1-x0)^2/(x2-2*x1+x0). (calculamos \(y_n=x_0^{(n)}\))

• • If |x-x0|<TOL then (miramos si el valor \(y_n\) cumple la condición de la tolerancia permitida)
    • Print(x); (damos el valor de la aproximación de cero o del punto fijo)
      
    • STOP. (Acabamos)
  • Set n++. (si no acabamos aumentar el contador del número de iteraciones en \(1\))
    
  • Set x0=x. (actualizamos el valor de \(x_0\))
    
• Print('El método ha fallado después de N iteraciones') (si hemos llegado hasta aquí significa que hemos alcanzado el número máximo de iteraciones sin que se cumpla la condición de la tolerancia)

Un polinomio de grado \(n\) tiene una expresión de la forma siguiente: \[ P_n(x)=a_n x^n+a_{n-1}x^{n-1}+\cdots +a_1 x+a_0, \] con \(a_n\neq 0\) para que tenga grado \(n\).

Los valores \(a_i\), \(i=0,1,\ldots,n\) se denominan coeficientes del polinomio \(P_n\).

El grado de \(P_n\) es el monomio (\(a_i x^i\)) de exponente \(i\) mayor con \(a_i\neq 0\). Si \(a_n\neq 0\), el grado de \(P_n\) será \(n\). Dicho \(a_n\) se llama coeficiente principal de \(P_n\).

Todos los métodos vistos hasta ahora se pueden aplicar al cálculo de ceros de polinomios pero necesitamos conocer un método eficiente para evaluar un polinomio y su derivada en un valor cualquiera \(z\) ya que, como veremos más adelante, si usamos su expresión teórica, el error cometido se dispara.

Primeramente, recordemos dos resultados de álgebra sobre los polinomios:

Teorema fundamental del álgebra. Sean \(P_n(x)\) un polinomio de grado \(n\) con coeficientes complejos. Entonces la ecuación \(P(x)=0\) siempre tiene una solución compleja.

La demostración del Teorema anterior se sale de los objetivos de curso actual de Métodos Numéricos.

Corolario. Sean \(P_n(x)\) un polinomio de grado \(n\) con coeficientes complejos y coeficiente principal \(a_n\). Entonces existen \(k\) valores complejos \(\alpha_1,\ldots,\alpha_k\) y \(k\) naturales \(m_1,\ldots,m_k\) tal que: \[ P_n(x)=a_n(x-\alpha_1)^{m_1}(x-\alpha_2)^{m_2}\cdots (x-\alpha_k)^{m_k}, \] con \(\displaystyle\sum_{i=1}^k m_i=n\).

Observación. Los valores \(\alpha_i\), \(i=1,\ldots,k\) son los ceros de multiplicidades \(m_1,\ldots,m_k\), respectivamente, del polinomio \(P_n(x)\).

Corolario. Sean \(P(x), Q(x)\) polinomios de grado como máximo \(n\) tal que existen \(k\) valores complejos con \(k>n\), \(x_1,\ldots, x_n\) tal que \(P(x_i)=Q(x_i)\) para \(i=1,\ldots,k\). Entonces \(P(x)=Q(x)\) para todo valor complejo \(x\).

Dado un polinomio \(P(x)\) de la forma \[ P_n(x)=a_n x^n+a_{n-1}x^{n-1}+\cdots +a_1 x+a_0, \] y dado un valor \(z\) real o complejo, para evaluar \(P(z)\) no debemos usar la expresión anterior ya que elevar un número \(x\) a una potencia cualquiera \(i\) hace aumentar el error cometido.

Veamos por qué. Consideramos la función \(g(x)=x^i\) y sea \(e_a(z)\) el error absoluto del número \(z\).

Aplicando la fórmula de propagación del error a la función \(g(x)\) obtenemos: \[ |g(z+e_a(z))-g(z)|\approx |g'(z)|\cdot e_a(z)=|i\cdot z^{i-1}|\cdot e_a(z). \] Es decir, el error al evaluar \(z^i\) es aproximadamente igual al error de \(z\) multiplicado por \(i\cdot |z|^{i-1}\). Si \(|z|>1\) y la potencia es grande, el error anterior se magnifica.

El método que debemos usar para evaluar \(P(z)\) se denomina regla de Horner y se basa en la división del polinomio \(P(x)\) por \(x-z\) y mirar el resto de la división: \[ P(x)=(x-z)\cdot Q(x)+R, \] donde \(R\) es el resto de división anterior.

Evaluando la expresión anterior para \(x=z\), obtenemos: \(P(z)=R\). En otras palabras, el valor de \(P(z)\) es precisamente el resto de la división anterior.

Para realizar la división de \(P(x)\) entre \(x-z\) podemos usar la regla de Ruffini:

El polinomio cociente \(Q(x)\) será: \[ Q(x)=b_{n-1}x^{n-1}+b_{n-2}x^{n-2}+\cdots + b_0, \] y el resto \(R\) de la división anterior será el valor de \(P(z)\).

Los coeficientes \(b_i\), \(i=0,\ldots, n-1\) de \(Q(x)\) se calculan de forma recurrente en función de los coeficientes \(a_i\), \(i=0,\ldots, n\) de \(P(x)\) de la manera siguiente: \[ b_{n-1}=a_n,\ b_i=a_{i+1}+b_{i+1}z,\ i=n-2,\ldots,0. \]

De esta forma, evitamos elevar un número a una potencia.

El método de Horner no sólo nos permite calcular \(P(z)\) sino también \(P'(z)\). ¿Cómo calculamos \(P'(z)\)?

Pues aplicando el método de Horner una vez más al polinomio cociente \(Q(x)\).

Es decir que si dividimos \(Q(x)\) por \(x-z\), el resto de dicha división es precisamente \(P'(z)\). Veamos por qué.

Recordemos que \[ P(x)=(x-z)\cdot Q(x)+P(z). \] Si derivamos la expresión anterior (respecto \(x\)), obtenemos: \[ P'(x)=Q(x)+(x-z)\cdot Q'(x), \] y evaluando en \(x=z\), obtenemos \(P'(z)=Q(z)\), donde \(Q(z)\) sería el valor del resto al dividir \(Q(x)\) por \(x-z\).

• INPUT n, a0,...,an,z. (Damos el grado del polinomio, sus coeficientes y el valor \(z\) en el que queremos calcular \(P_n(z)\))
  
• Set b=an. (En \(b\) iremos almacenando los coeficientes \(b_n\) de la división de \(P(x)\) entre \(x-z\) y en \(c\), los coeficientes de las divisiones de \(Q(x)\) entre \(x-z\))
• Set c=an.
• For j=n-1,...,1 (Calculamos los demás coeficientes)
  • Set b=z*b+aj.
    
  • Set c=z*c+b.
    
• Set b=z*b+a0. (Calculamos \(P(z)\))
• Print b,c. (Damos los valores de \(P(z)\) y \(P'(z)\) respectivamente)

Observación. Podemos aplicar el método de Horner que nos da el valor de \(P(z)\) y \(P'(z)\) para cualquier valor de \(z\), al método de Newton-Raphson para el cálculo de ceros de polinomios.

Una vez hallado un cero \(\hat{x}\) del polinomio \(P_n(x)\) de grado \(n\), podemos considerar el polinomio \(P_{n-1}=\frac{P_n(x)}{x-\hat{x}}\) de grado \(n-1\) que tiene los mismos ceros que \(P_n(x)\) excepto el cero hallado \(\hat{x}\).

Entonces podemos intentar hallar algun cero de \(P_{n-1}(x)\) y realizar el mismo proceso anterior.

De esta forma, podemos llegar a calcular todos los ceros de \(P_n(x)\).

A dicho proceso se le denomina deflación.

El problema que nos planteamos en esta sección es el cálculo de ceros complejos para polinomios \(P_n(x)\) de grado \(n\) con coeficientes reales.

Si aplicamos los métodos vistos anteriormente, bisección, Newton-Raphson, secante o regula-falsi y empezamos con un valor inicial \(x_0\) real, los valores siguientes de la sucesión serán todos reales y será imposible hallar un cero complejo.

Una solución es empezar con un valor complejo, de esta forma los valores obtenidos serían complejos y podríamos llegar al cero.

Uno de los motivos del comportamiento que observamos en el ejemplo anterior nos lo da la proposición siguiente:

Proposición. Sea \(P_n(x)=a_n x^n+a_{n-1}x^{n-1}+\cdots +a_1 x+a_0\) un polinomio con coeficientes reales, \(a_i\in\mathbb{R}\), \(i=0,1,\ldots,n\). Sea \(z=\alpha +\beta\mathrm{i}\) una raíz compleja de \(P_n(x)\). Entonces el conjugado de \(z\), \(\overline{z}=\alpha -\beta\mathrm{i}\) también es un cero del polinomio \(P_n(x)\).

Entonces se hace necesario tener un método para hallar ceros complejos de polinomios con una velocidad de convergencia lo más rápida posible.

El método de Müller es un buen candidato para realizar la tarea anterior.

El método de Müller usa una estrategia parecida al método de la secante.

Recordemos que en el método de la secante, para hallar el cero de la función \(f\) o resolver \(f(x)=0\), calculábamos la recta secante que pasaba por los puntos \((x_{n-2},f(x_{n-2}))\) y \((x_{n-1},f(x_{n-1}))\) y hallábamos el corte de dicha recta con el eje \(X\) para hallar el punto siguiente de la sucesión \(x_n\).

En el método de Müller, hacemos una cosa parecida, pero en lugar de consirar la recta secante anterior, consideramos la parábola que pasa por los puntos \((x_{n-3},f(x_{n-3}))\), \((x_{n-2},f(x_{n-2}))\) y \((x_{n-1},f(x_{n-1}))\) y hallamos el corte de dicha parábola con el eje \(X\) para hallar el siguiente punto de la sucesión \(x_n\).

Observar que los valores \(x_i\) y \(f(x_i)\) son complejos.

Más concretamente, sea \(P_2(x)\) la parábola anterior que escribimos de la forma siguiente: \[ P_2(x)=a (x-x_{n-1})^2+b (x-x_{n-1})+c. \] Para hallar los coeficientes \(a\), \(b\) y \(c\) hemos de imponer que la parábola anterior pasa por los puntos \((x_{n-3},f(x_{n-3}))\), \((x_{n-2},f(x_{n-2}))\) y \((x_{n-1},f(x_{n-1}))\): \[ \left. \begin{align*} f(x_{n-3}) = & a (x_{n-3}-x_{n-1})^2+b (x_{n-3}-x_{n-1})+c, \\ f(x_{n-2}) = & a (x_{n-2}-x_{n-1})^2+b (x_{n-2}-x_{n-1})+c, \\ f(x_{n-1}) = & a (x_{n-1}-x_{n-1})^2+b (x_{n-1}-x_{n-1})+c, \end{align*} \right\} \]

Los valores \(a\), \(b\) y \(c\) se calculan resolviendo el sistema anterior: \[ {\scriptsize \begin{align*} c= & f(x_{n-1}),\\ b = & \frac{(x_{n-3}-x_{n-1})^2 (f(x_{n-2})-f(x_{n-1}))-(x_{n-2}-x_{n-1})^2 (f(x_{n-3})-f(x_{n-1}))}{(x_{n-3}-x_{n-1})(x_{n-2}-x_{n-1})(x_{n-3}-x_{n-2})},\\ a = & \frac{(x_{n-2}-x_{n-1}) (f(x_{n-3})-f(x_{n-1}))-(x_{n-3}-x_{n-1}) (f(x_{n-2})-f(x_{n-1}))}{(x_{n-3}-x_{n-1})(x_{n-2}-x_{n-1})(x_{n-3}-x_{n-2})}. \end{align*} } \] Para hallar \(x_n\) debemos resolver \(P_2(x)=a (x-x_{n-1})^2+b (x-x_{n-1})+c=0\) cuya solución es: \[ x-x_{n-1}=\frac{-b\pm \sqrt{b^2-4ac}}{2a},\ \Rightarrow x_n=x=x_{n-1}+\frac{-b\pm \sqrt{b^2-4ac}}{2a}. \]

Vemos que tenemos que elegir entre \(\pm\) en la expresión de \(x_n\). Elegimos el signo que hace que \(x_n\) esté más cerca de \(x_{n-1}\).

Es decir,

• si \(\left|\frac{-b+ \sqrt{b^2-4ac}}{2a}\right|<\left|\frac{-b - \sqrt{b^2-4ac}}{2a}\right|\), elegimos el signo \(+\),
• en caso contrario, elegimos el signo \(-\).

• INPUT x0, x1, x2, TOL, N. (Damos los valores de \(x_0\), \(x_1\), \(x_2\), la tolerancia permitida TOL y el número máximo de iteraciones \(N\))
• Set h1=x1-x0. (Introducimos los \(h_i\) como las diferencias entre dos valores iniciales consecutivos)
• Set h2=x2-x1.
• Set y1=(f(x1)-f(x0))/h1. (Introducimos las \(y_i\) que serán los cocientes incrementales que aparecen en las expresiones de \(a\), \(b\) y \(c\))
• Set y2=(f(x2)-f(x1))/h2.
• Set a=(y2-y1)/(h2+h1). (Calculamos el valor de \(a\))
• Set n=3 (n será un contador que nos dirá el valor de \(n\) en la sucesión \(x_n\))

• While n <= N do (Mientras no llegemos al máximo de iteraciones hacer lo siguiente)
  • Set b=y2+h2*a. (Calculamos el valor de \(b\))
    
  • Set D=sqrt(b^2-4*f(x2)*a). (Calculamos la raíz del discriminante o \(\sqrt{b^2-4ac}\))
    
  • If |-b+D|< |-b-D| then (hemos de elegir el signo \(+\) en la expresión de \(x_n\))
    • Set E=-b+D.
      
  • Else (hemos de elegir el signo \(-\) en la expresión de \(x_n\))
    • Set E=-b-D.
  • Set x=x2+E/(2*a). (Calculamos el valor de \(x_n\))

• • If |x-x2|<TOL then (Miramos si el valor \(x_n\) cumple la condición de la tolerancia permitida)
    • Print (x). (Damos \(x\) como el valor aproximado del cero buscado)
      
    • STOP. (Paramos)
  • Set x0=x1. (Actualizamos los valores para la próxima iteración)
    
  • Set x1=x2.
  • Set x2=x.
  • Set h1=x1-x0.
  • Set h2=x2-x1.

• • Set y1=(f(x1)-f(x0))/h1.
    
  • Set y2=(f(x2)-f(x1))/h2.
  • Set a=(y2-y1)/(h2+h1).
  • Set n++.
• Print('El método ha fallado después de N iteraciones') (si hemos llegado hasta aquí significa que hemos alcanzado el número máximo de iteraciones sin que se cumpla la condición de la tolerancia)