Orden de convergencia
Convergencia lineal, cuadrática y superlineal, y cómo leer el orden en la sucesión de errores.
Recordemos el análisis de la serie de errores del método de punto fijo, que vimos en la nota de raíces y puntos fijos, y usémoslo para darle nombre a la velocidad con que converge un método. Después lo aplicamos a Newton, que vimos en la nota anterior, y a la secante. Seguimos las secciones 4.2, 4.3 y 4.4 del libro de Driscoll y Braun.
1 El error en la iteración de punto fijo
Sea \(p\) la solución exacta, es decir, \(g(p)=p\) o, de forma equivalente, \(f(p)=g(p)-p=0\). El error absoluto en el paso \(k\) es \(\varepsilon_k=x_k-p\). Como \(x_{k+1}=g(x_k)\), expandiendo \(g\) en Taylor alrededor de \(p\) y usando \(g(p)=p\),
\[\varepsilon_{k+1}+p=g(p+\varepsilon_k)=g(p)+g'(p)\,\varepsilon_k+o(\varepsilon_k)=p+g'(p)\,\varepsilon_k+o(\varepsilon_k).\]
Para \(\varepsilon_k\) pequeño el resto \(o(\varepsilon_k)\) se puede despreciar, y queda
\[\boxed{\varepsilon_{k+1}\approx g'(p)\,\varepsilon_k.}\]
Ya vimos que para que el método funcione cerca de una semilla se requiere \(|g'(p)|<1\), y que entonces el error decrece como un múltiplo escalar del error anterior. Esto motiva la definición.
2 Convergencia lineal
Si \(x_k\to x^*\) y el error \(\varepsilon_k=x_k-x^*\) satisface
\[\lim_{k\to\infty}\frac{|\varepsilon_{k+1}|}{|\varepsilon_k|}=\sigma<1,\]
se dice que la convergencia es lineal.
En el punto fijo, \(\sigma=|g'(p)|\). Supongamos que la razón es exactamente \(\sigma\) en cada paso. Entonces
\[|\varepsilon_k|=\sigma\,|\varepsilon_{k-1}|=\sigma^2\,|\varepsilon_{k-2}|=\cdots=\sigma^k\,|\varepsilon_0|=\sigma^k c,\qquad c=|\varepsilon_0|,\]
de donde, al tomar logaritmos,
\[\log|\varepsilon_k|=k\log\sigma+\log c.\]
Es una relación lineal entre \(\log|\varepsilon_k|\) y \(k\), con pendiente \(\log\sigma<0\): si graficamos los errores en escala logarítmica contra \(k\), aparece una recta. De ahí el nombre.
Por ejemplo, para \(g(x)=\cos x\) y \(x_0=1\), el punto fijo es \(p\approx0{,}7390851\) y \(g'(p)=-\operatorname{sen}(p)\approx-0{,}67361\). Los cocientes \(\varepsilon_{k+1}/\varepsilon_k\) se acercan a ese valor, alternando de signo porque \(g'(p)<0\), de modo que \(\sigma\approx0{,}6736\):
| \(k\) | \(\varepsilon_k\) | \(\varepsilon_{k+1}/\varepsilon_k\) |
|---|---|---|
| \(0\) | \(0{,}26091\) | \(-0{,}76187\) |
| \(5\) | \(-0{,}037716\) | \(-0{,}65952\) |
| \(10\) | \(0{,}0051522\) | \(-0{,}67551\) |
| \(20\) | \(9{,}9267\times10^{-5}\) | \(-0{,}67365\) |
| \(30\) | \(1{,}9095\times10^{-6}\) | \(-0{,}67361\) |
3 Convergencia cuadrática: el método de Newton
Veamos ahora la convergencia del método de Newton con este lenguaje del error; en la nota anterior vimos que era cuadrática. Sea \(r\) una raíz de \(f\) con \(f'(r)\ne0\), y \(f\) dos veces diferenciable cerca de \(r\). Con \(\varepsilon_k=x_k-r\), restando \(r\) a ambos lados de \(x_{k+1}=x_k-f(x_k)/f'(x_k)\) se obtiene
\[\varepsilon_{k+1}=\varepsilon_k-\frac{f(r+\varepsilon_k)}{f'(r+\varepsilon_k)}.\]
Hacemos una expansión de Taylor de orden 2 en \(f\) y de orden 1 en \(f'\). Como \(f(r)=0\),
\[\varepsilon_{k+1}=\varepsilon_k-\frac{f'(r)\,\varepsilon_k+\tfrac12f''(r)\,\varepsilon_k^2+o(\varepsilon_k^2)}{f'(r)+f''(r)\,\varepsilon_k+o(\varepsilon_k)}=\varepsilon_k-\varepsilon_k\,\frac{1+\dfrac{f''(r)}{2f'(r)}\,\varepsilon_k+o(\varepsilon_k)}{1+\dfrac{f''(r)}{f'(r)}\,\varepsilon_k+o(\varepsilon_k)}.\]
El denominador es \(1+\dfrac{f''(r)}{f'(r)}\varepsilon_k+o(\varepsilon_k)\), con \(\left|\dfrac{f''(r)}{f'(r)}\varepsilon_k\right|<1\) porque \(\varepsilon_k\) es pequeño. Por la serie geométrica,
\[\frac{1}{1+\dfrac{f''(r)}{f'(r)}\,\varepsilon_k+o(\varepsilon_k)}=1-\frac{f''(r)}{f'(r)}\,\varepsilon_k+o(\varepsilon_k),\]
de donde
\[\begin{aligned} \varepsilon_{k+1}&=\varepsilon_k-\varepsilon_k\left(1+\frac{f''(r)}{2f'(r)}\,\varepsilon_k+o(\varepsilon_k)\right)\left(1-\frac{f''(r)}{f'(r)}\,\varepsilon_k+o(\varepsilon_k)\right)\\ &=\varepsilon_k-\varepsilon_k\left(1-\frac{f''(r)}{2f'(r)}\,\varepsilon_k+o(\varepsilon_k)\right)\\ &=\frac12\,\frac{f''(r)}{f'(r)}\,\varepsilon_k^2+o(\varepsilon_k^2). \end{aligned}\]
Si \(x_k\to x^*\) y el error \(\varepsilon_k=x_k-x^*\) satisface
\[\lim_{k\to\infty}\frac{|\varepsilon_{k+1}|}{|\varepsilon_k|^2}=L>0,\]
se dice que la convergencia es cuadrática.
Es el caso de Newton: dividiendo el resultado anterior por \(\varepsilon_k^2\),
\[\frac{|\varepsilon_{k+1}|}{|\varepsilon_k|^2}=\frac12\left|\frac{f''(r)}{f'(r)}\right|+o(1),\]
y tomando \(k\to\infty\) a ambos lados, \(L=\dfrac12\dfrac{|f''(r)|}{|f'(r)|}\), que es positivo siempre que \(f''(r)\ne0\). Tomando logaritmos en \(|\varepsilon_{k+1}|\approx L\,|\varepsilon_k|^2\),
\[\log|\varepsilon_{k+1}|\approx2\log|\varepsilon_k|+\log L,\qquad\frac{\log|\varepsilon_{k+1}|}{\log|\varepsilon_k|}\approx2+\frac{\log L}{\log|\varepsilon_k|}\longrightarrow2\quad(k\to\infty),\]
porque \(\log|\varepsilon_k|\to-\infty\). En palabras, cada paso de Newton duplica aproximadamente el número de cifras correctas.
Si \(x_k\to x^*\) y existen constantes \(\alpha>1\) y \(L>0\) tales que
\[\lim_{k\to\infty}\frac{|\varepsilon_{k+1}|}{|\varepsilon_k|^\alpha}=L,\]
se dice que la convergencia es superlineal de orden \(\alpha\). Con \(\alpha=2\) es cuadrática, y con \(\alpha>2\) se llama supercuadrática.
4 La secante: superlineal, pero no cuadrática
En el método de Newton reemplacemos \(f'(x_k)\) por una pendiente \(m\) que aproxime la tangente: la de la recta secante por los dos últimos puntos de la gráfica,
\[m=\frac{f(x_k)-f(x_{k-1})}{x_k-x_{k-1}}.\]
El método de la secante es entonces
\[x_{k+1}=x_k-\frac{f(x_k)}{m}=x_k-f(x_k)\,\frac{x_k-x_{k-1}}{f(x_k)-f(x_{k-1})},\]
y necesita dos iterados iniciales. Geométricamente, \(x_{k+1}\) es el punto donde la recta secante corta el eje horizontal.
4.1 El error de la secante
Con \(\varepsilon_j=x_j-r\) se tiene \(x_k-x_{k-1}=\varepsilon_k-\varepsilon_{k-1}\), y restando \(r\) en la fórmula anterior,
\[\varepsilon_{k+1}=\varepsilon_k-f(x_k)\,\frac{\varepsilon_k-\varepsilon_{k-1}}{f(x_k)-f(x_{k-1})}.\]
Llevando todo al denominador común \(f(x_k)-f(x_{k-1})\), los términos \(\varepsilon_k f(x_k)\) se cancelan:
\[\varepsilon_{k+1}=\frac{\varepsilon_k\,[f(x_k)-f(x_{k-1})]-f(x_k)\,(\varepsilon_k-\varepsilon_{k-1})}{f(x_k)-f(x_{k-1})}=\frac{f(x_k)\,\varepsilon_{k-1}-f(x_{k-1})\,\varepsilon_k}{f(x_k)-f(x_{k-1})}.\]
Expandimos ahora con Taylor de orden 2 alrededor de \(r\): \(f(x_j)=f'(r)\,\varepsilon_j+\tfrac12f''(r)\,\varepsilon_j^2+o(\varepsilon_j^2)\) para \(j=k,k-1\). En el numerador, los términos con \(f'(r)\) se cancelan entre sí y queda
\[f(x_k)\,\varepsilon_{k-1}-f(x_{k-1})\,\varepsilon_k=\tfrac12f''(r)\,\varepsilon_k^2\varepsilon_{k-1}-\tfrac12f''(r)\,\varepsilon_{k-1}^2\varepsilon_k+\cdots=\tfrac12f''(r)\,\varepsilon_k\varepsilon_{k-1}\,(\varepsilon_k-\varepsilon_{k-1})+\cdots,\]
mientras que el denominador es \(f'(r)\,(\varepsilon_k-\varepsilon_{k-1})+\cdots\). Por lo tanto
\[\varepsilon_{k+1}\approx\frac{\tfrac12f''(r)\,\varepsilon_k\varepsilon_{k-1}\,(\varepsilon_k-\varepsilon_{k-1})}{f'(r)\,(\varepsilon_k-\varepsilon_{k-1})}=\frac12\,\frac{f''(r)}{f'(r)}\,\varepsilon_k\varepsilon_{k-1}.\]
Los puntos suspensivos también llevan el factor \(\varepsilon_k-\varepsilon_{k-1}\): un término \(c_i\varepsilon^i\) del desarrollo de \(f\), con \(i\ge2\), aporta \(c_i\,\varepsilon_k\varepsilon_{k-1}\,(\varepsilon_k^{\,i-1}-\varepsilon_{k-1}^{\,i-1})\) al numerador, y con \(i\ge1\) aporta \(c_i\,(\varepsilon_k^{\,i}-\varepsilon_{k-1}^{\,i})\) al denominador. Ambos son divisibles por \(\varepsilon_k-\varepsilon_{k-1}\), y por eso la cancelación es legítima.
\[\boxed{\varepsilon_{k+1}\approx\frac12\,\frac{f''(r)}{f'(r)}\,\varepsilon_k\varepsilon_{k-1}}\]
Si \(\varepsilon_{k-1}\) fuera igual a \(\varepsilon_k\), la secante sería la tangente y recuperaríamos la fórmula de Newton, \(\varepsilon_{k+1}\approx\frac12\frac{f''(r)}{f'(r)}\varepsilon_k^2\).
4.2 El orden de la secante
Como \(\varepsilon_{k+1}\) depende de los dos errores anteriores, no sabemos de antemano a qué potencia de \(\varepsilon_k\) es proporcional. Supongamos que los errores tienen esta forma, con \(c>0\) y \(\alpha>0\):
\[|\varepsilon_{k+1}|=c\,|\varepsilon_k|^\alpha,\qquad|\varepsilon_k|=c\,|\varepsilon_{k-1}|^\alpha.\]
Sea \(C=\left|\dfrac{f''(r)}{2f'(r)}\right|\). Por la fórmula del error de la secante,
\[|\varepsilon_{k+1}|\approx C\,|\varepsilon_k|\,|\varepsilon_{k-1}|=C\,c\,|\varepsilon_{k-1}|^{\alpha+1},\]
y por otro lado
\[|\varepsilon_{k+1}|=c\,|\varepsilon_k|^\alpha=c\,\big(c\,|\varepsilon_{k-1}|^\alpha\big)^\alpha=c^{1+\alpha}\,|\varepsilon_{k-1}|^{\alpha^2}.\]
Igualando ambas expresiones,
\[|\varepsilon_{k-1}|^{\alpha^2-\alpha-1}\approx C\,c^{-\alpha}.\]
El lado derecho es una constante positiva, mientras que \(|\varepsilon_{k-1}|\to0\). Esto solo es posible si el exponente es cero, es decir, si
\[\alpha^2=\alpha+1\qquad\Longrightarrow\qquad\alpha=\frac{1+\sqrt5}{2}=\phi\approx1{,}618,\]
la razón áurea, que es la única raíz positiva. De paso, la constante queda determinada por \(c^\alpha=C\), o sea \(c=\left|\dfrac{f''(r)}{2f'(r)}\right|^{\alpha-1}\), porque \(1/\phi=\phi-1\).
Entonces los errores de la secante satisfacen
\[\lim_{k\to\infty}\frac{|\varepsilon_{k+1}|}{|\varepsilon_k|^\alpha}=L>0,\qquad\alpha=\frac{1+\sqrt5}{2}\in(1,2),\]
es decir, la convergencia es superlineal pero no cuadrática. Tomando logaritmos a ambos lados, \(\log|\varepsilon_{k+1}|\approx\alpha\log|\varepsilon_k|+\log L\), de donde
\[\frac{\log|\varepsilon_{k+1}|}{\log|\varepsilon_k|}\longrightarrow\alpha.\]
5 Leer el orden en la práctica
La sucesión de cocientes \(\log|\varepsilon_{k+1}|/\log|\varepsilon_k|\) sirve en general: su límite \(\alpha\) es el orden de convergencia del algoritmo que estemos probando. Para calcularla necesitamos errores, así que usamos aritmética de precisión extendida (200 dígitos) y una raíz de referencia calculada con mucha más precisión que los iterados. En doble precisión el error llega a \(10^{-16}\) en pocos pasos y la sucesión se estanca antes de mostrar la tendencia.
Tomamos \(f(x)=xe^x-2\), cuya raíz es \(r\approx0{,}8526055020137255\). Con Newton desde \(x_0=1\):
| \(k\) | \(\lvert\varepsilon_k\rvert\) | \(\log\lvert\varepsilon_{k+1}\rvert\big/\log\lvert\varepsilon_k\rvert\) | \(\lvert\varepsilon_{k+1}\rvert\big/\lvert\varepsilon_k\rvert^2\) |
|---|---|---|---|
| \(0\) | \(0{,}1474\) | \(2{,}184\) | \(0{,}70305\) |
| \(1\) | \(0{,}01527\) | \(2{,}065\) | \(0{,}76244\) |
| \(2\) | \(1{,}779\times10^{-4}\) | \(2{,}030\) | \(0{,}76980\) |
| \(3\) | \(2{,}436\times10^{-8}\) | \(2{,}015\) | \(0{,}76989\) |
| \(4\) | \(4{,}567\times10^{-16}\) | \(2{,}007\) | \(0{,}76989\) |
| \(5\) | \(1{,}606\times10^{-31}\) | \(2{,}004\) | \(0{,}76989\) |
El cociente de logaritmos tiende a \(2\), y el cociente \(|\varepsilon_{k+1}|/|\varepsilon_k|^2\) se estabiliza en \(0{,}76989\), que es justamente \(\frac12|f''(r)/f'(r)|\). Con la secante, partiendo de \(x_1=1\) y \(x_2=0{,}5\):
| \(k\) | \(\lvert\varepsilon_k\rvert\) | \(\log\lvert\varepsilon_{k+1}\rvert\big/\log\lvert\varepsilon_k\rvert\) |
|---|---|---|
| \(4\) | \(0{,}01303\) | \(1{,}787\) |
| \(5\) | \(4{,}275\times10^{-4}\) | \(1{,}594\) |
| \(6\) | \(4{,}270\times10^{-6}\) | \(1{,}649\) |
| \(7\) | \(1{,}405\times10^{-9}\) | \(1{,}619\) |
| \(8\) | \(4{,}620\times10^{-15}\) | \(1{,}625\) |
| \(9\) | \(4{,}999\times10^{-24}\) | \(1{,}620\) |
| \(10\) | \(1{,}778\times10^{-38}\) | \(1{,}620\) |
| \(11\) | \(6{,}845\times10^{-62}\) | \(1{,}619\) |
El cociente de logaritmos se acerca a \(\phi\approx1{,}618\). El cociente \(|\varepsilon_{k+1}|/|\varepsilon_k|^\alpha\) oscila alrededor de \(0{,}85\) y se acerca lentamente al valor que predice la deducción, \(\left|\frac{f''(r)}{2f'(r)}\right|^{\alpha-1}\approx0{,}8508\).
Para convergencia lineal el cociente de logaritmos también tiende a \(\alpha=1\), pero muy despacio: con \(|\varepsilon_k|=\sigma^kc\) vale \(1+\dfrac{\log\sigma}{k\log\sigma+\log c}\), que decrece como \(1/k\). En la iteración \(x_{k+1}=\cos x_k\) de arriba vale \(1{,}20\) en \(k=0\) y todavía \(1{,}02\) en \(k=39\). Para detectar una convergencia lineal conviene mirar directamente el cociente \(\varepsilon_{k+1}/\varepsilon_k\).
Resuelva todos los ejercicios de las secciones 4.2, 4.3 y 4.4 de Fundamentals of Numerical Computation (iteración de punto fijo, método de Newton y métodos basados en interpolación), disponibles en https://fncbook.com/julia/fixed-point/#exercises, https://fncbook.com/julia/newton/#exercises y https://fncbook.com/julia/secant/#exercises.