A continuación se presenta el código en Matlab® empleado para generar la figura anterior, excepto por los títulos x(k) que se añadieron con un editor de imágenes.
% macro para construir la Figura 2.4
f1=@(x) x-x^2+0.5*x^3; % función no lineal f1(x)
f1p=@(x) 1-2*x+1.5*x^2; % derivada de la función f1’(x)
fplot(f1,[-2 4],’k-’);
axis square;
xlabel(‘x’);
ylabel(‘y’);
hold on;
plot([-2 4],[0 0],’k-’)
% secuencia iterativa del método de Newton
x(1)=3.5;
f(1)=f1(x(1));
for i=2:4,
% iteración del método de Newton
x(i)=x(i-1)-f(i-1)/f1p(x(i-1));
f(i)=f1(x(i));
% secuencia de gráficos del método de Newton
xpl=[x(i-1) x(i-1)];
fpl=[0 f(i-1)];
plot(xpl,fpl,’k:’);
xpl=[x(i-1) x(i)];
fpl=[f(i-1) 0];
plot(xpl,fpl,’k:’);
end
Notas:
i) La función de actualización es: φ(x) = x - f(x)/f’(x)
ii) El método es de segundo orden: cuadráticamente convergente. Esto significa que el número de cifras decimales correctas se duplica en cada iteración cuando x(k) → x*
Podemos apreciar en la ecuación 2.7 que define la actualización del método de Newton, que se pueden presentar problemas si cerca de la solución buscada la derivada f’ se hace cero o muy pequeña. El siguiente teorema nos aclara qué se puede esperar de la iteración 2.7.
Teorema: Si la solución x* de f(x) = 0 está en un intervalo I en donde f’(x) ≠ 0 y además f’ es continua en I y se cumple la condición de convexidad:
Entonces para todo x(0) ∈ I, la iteración de Newton produce una secuencia {x(k)} que satisface exactamente una de las siguientes condiciones:
i) x(k) ∈ I y la secuencia {x(k)} converge monótonamente a x*
ii) f(x(0)) f(x(1)) < 0 y la secuencia {x(k)} converge monótonamente a x*
iii) x(1)
Corolario 1: La condición de convexidad se puede cambiar a una de concavidad:
Y el teorema sigue siendo válido.
Corolario 2: En vez de condiciones de convexidad o concavidad, se puede utilizar la condición: f” (x) ≠ 0 y es continua ∀ x ∈ I.
2.3.2 Interpolación cuadrática
Dados los problemas del método de Newton clásico cuando la derivada cambia de signo y la iteración se torna inestable, se han propuesto esquemas para evitar dividir por la derivada f’(x).
En la metodología de interpolación cuadrática se aproxima f(x) en torno a x = x(k) con un polinomio de segundo orden:
y x(k+1) se obtiene haciendo P(x(k+1)) = 0:
Puesto que se supone que x(k+1) ≈ x(k), se escoge la raíz con menor valor absoluto en la fórmula de iteración.
Comentarios:
i) La cancelación de términos de similar valor absoluto, pero distinto signo, puede conducir a pérdida de cifras significativas en la aritmética.
ii) Dado que no se puede garantizar en general que el discriminante γ(k) en la ecuación 2.11 sea positivo, lo que se hace es remplazar uno de los factores (x-x(k)) por la expresión del método de Newton clásico 2.7, con lo que se obtiene el método de Newton extendido:
El siguiente ejemplo muestra que este método permite reducir las oscilaciones que ocurren con el método de Newton clásico, aunque no se puede concluir que el algoritmo 2.12 sea siempre robusto con respecto al estimador inicial x(1).
Ejemplo 2.5. Comparación del método de Newton clásico y su versión extendida
Compare el desempeño del método de Newton clásico con su variante extendida para el caso de f(x) = x3-x-1 partiendo con x(1) = -1.
Se obtiene la siguiente figura, donde se aprecia que el método clásico oscila en las primeras 10 iteraciones sin acercarse a la raíz correcta, mientras que al método extendido le toma solamente 4 o 5 iteraciones llegar a la raíz correcta, aun cuando la derivada cerca de x(1) = -1 va cambiando de signo.
FIGURA 2.5. Desempeño de la iteración de Newton clásica y extendida
A continuación se presenta el código de Matlab® que permite hacer los cálculos y la figura anterior.
% macro que resuelve el ejemplo 2.4
f=@(x) x^3-x-1;
fp=@(x) 3*x^2-1;
fpp=@(x) 6*x;
% secuencia iterativa del método de Newton y de su variante cuadrática
xN=zeros(10,1);
fN=zeros(10,1);
xN(1)=-1;
fN(1)=f(xN(1));
xC=zeros(10,1);
fC=zeros(10,1);
xC(1)=-1;
fC(1)=f(xC(1));
for i=2:10,
% iteración del método de Newton
xN(i)=xN(i-1)-fN(i-1)/fp(xN(i-1));
fN(i)=f(xN(i));
% iteración por método de Newton extendido
der=fp(xC(i-1));
xC(i)=xC(i-1)-fC(i-1)/(der-0.5*fC(i-1)*fpp(xC(i-1))/der);
fC(i)=f(xC(i));
end
% gráfico de comparación de resultados
plot(1:10,xN,’k-o’,1:10,xC,’k--^’);
axis square;
title(‘f(x) = x^3 - x - 1’);
xlabel(‘Número de iteraciones’);
ylabel(‘Raíz estimada’);
legend(‘Newton clásico’,’Newton extendido’);
dc(3)