Método de Newton-Raphson y algunas propiedades

El método de Newton-Raphson es un método abierto para resolver ecuaciones no lineales. Resolver se entiende como encontrar los valores de x donde la función f(x) es cero

El valor particular de x donde la función es cero se conoce como raíz.
El método de Newton-Raphson se deriva de la serie de Taylor, donde se aproxima una función por medio de la serie y esa serie se iguala a cero, para encontrar la aproximación a la raíz se despeja de la serie truncada de primer orden

Donde x0 es un valor inicial con la que se hace la primer aproximación y f'(x) es la derivada de f(x).
El algoritmo es sencillo:

function [x_r,k] = nr(datos)
f   = datos.funcion;
df  = datos.derivada;
x   = datos.x0;
n   = datos.niter;
tol = datos.tol;
k   = 0;
x_r = NaN;
while k<n
    x = x - f(x)/df(x);
    k = k + 1;
    if abs(f(x))<tol
        x_r = x;
        break;
    end
end
end

Donde iteraciones es el número máximo de iteraciones en el ciclo, k es el contador de iteraciones, tolerancia es el valor máximo que puede tener la función para aceptarse como la raíz de f(x). Ambos valores son arbitrarios y se dan de acuerdo al problema que se quiere resolver.
Por ejemplo, si se quiere saber la raíz de la siguiente función

obtenemos su derivada que es

Podemos dar un valor inicial por ejemplo x0=1, un valor máximo de iteraciones = 100 y una tolerancia = 1e-5 . Todo eso lo creamos en una estructura

datos.funcion  = @(x) x.^2-2;
datos.derivada = @(x) 2*x;
datos.x0       = 1;
datos.niter    = 100;
datos.tol      = 1e-5;

Ejecutamos el programa de Newton-Raphson para encontrar la raíz

[raiz,k]=nr(datos)
raiz = 1.4142
k = 3

El resultado es raiz=1.4142 en 3 iteraciones. Si cambiamos el valor inicial por ejemplo x0=0.5 obtenemos

datos.x0       = 0.5;
[raiz,k]=nr(datos)
raiz = 1.4142
k = 5

Observamos que llegamos al mismo resultado de la raíz pero con 5 iteraciones, ahora probemos con otro valor de x0=0.1

datos.x0       = 0.1;
[raiz,k]=nr(datos)
raiz = 1.4142
k = 7

Sigue obteniendo el resultado pero ahora en 7 iteraciones, hagamos otra prueba, ahora con x0=0.001

datos.x0       = 0.001;
[raiz,k]=nr(datos)
raiz = 1.4142
k = 14

Observamos que igualmente llega al resultado pero ahora le toma 14 iteraciones, notamos que si el valor inicial se acerca a 0 entonces le toma más iteraciones llegar a la raíz, no es posible que x0=0 porque la derivada es cero y el valor de x se indetermina y no es posible obtener la raíz.

datos.x0       = 0;
[raiz,k]=nr(datos)
raiz = NaN
k = 100

Podemos pensar que existen valores de x0 donde es más fácil encontrar la raíz y existen otros donde encontrar la raíz es un poco más tardado.
Seamos un poco curiosos y aberigüemos qué comportamiento tiene el método de Newton-Raphson para otros valores iniciales. Si hacemos un proceso donde cambiemos el valor inicial de x0 en el intervalo [-2, 2] y obtenemos el número de iteraciones k que le toma llegar a la raíz en cada caso y grafiquemos esos resultados para observar el comportamiento

n=100;
nPuntos=n;
rangox=[-2 2];
xs=linspace(rangox(1),rangox(2),n);
ks=zeros(n,1);
for j=1:nPuntos
    datos.x0=xs(j);
    [xr,ks(j)]=nr(datos);
end
plot(xs,ks)
title('Newton-Raphson para f(x)=x^2-2')
xlabel('x inicial')
ylabel('iteraciones')
grid

Observamos dos mínimos cercanos a -1.4142 y 1.4142 , esto se debe a que existen dos raíces para esta función.
Será el mismo comportamiento para otras funciones?
Cambiemos ahora la función como

datos.funcion  = @(x) x.^3-1;
datos.derivada = @(x) 3*x.^2;

y hagamos el mismo ejercicio para observar el comportamiento

n=1000;
nPuntos=n;
rangox=[-2 2];
xs=linspace(rangox(1),rangox(2),n);
ks=zeros(n,1);
for j=1:nPuntos
    datos.x0=xs(j);
    [xr,ks(j)]=nr(datos);
end
plot(xs,ks)
title('Newton-Raphson para f(x)=x^3-1')
xlabel('x inicial')
ylabel('iteraciones')
grid

A diferencia de la función cuadrática, ésta función al ser cúbica tiene 3 raíces, una es real y las otras complejas. El método de Newton-Raphson es capaz de encontrar raíces complejas, solo que el valor inicial x0 debe ser un valor complejo. Entonces debemos crear un valor complejo para los valores iniciales, y debemos variar tanto el valor real como el imaginario de cada valor complejo y hacer todas las combinaciones

xs=linspace(rangox(1),rangox(2),n); %valores reales
ys=linspace(rangoy(1),rangoy(2),n); %valores imaginarios
[x, y] = meshgrid(xs, ys); %combinaciones de cada valor

Para poder graficar tanto el valor real como el imaginario y además el número de iteraciones que requiere el método para encontrar la raíz tendremos que usar una gráfica en tres dimensiones donde el plano x-y es el valor complejo y el eje z son las iteraciones.

n=60;
nPuntos=n*n;
lienzo = zeros(n, n);
rangox=[-2 2];
rangoy=[-2 2];
xs=linspace(rangox(1),rangox(2),n); %valores reales
ys=linspace(rangoy(1),rangoy(2),n); %valores imaginarios
[x, y] = meshgrid(xs, ys); %combinaciones de cada valor
ks=zeros(n,1);
for j=1:nPuntos
    datos.x0=x(j)+i*y(j); %se crea el valor complejo con i
    [xr,ks(j)]=nr(datos);
    lienzo(j)=ks(j);
end
surfc(x,y,lienzo)
title('Newton-Raphon f(x)=x^3-1')
xlabel('Parte real')
ylabel('Parte imaginaria')
zlabel('Iteraciones')
view([-23.5 62.0])

Observamos zonas de convergencia que son las más bajas y zonas de divergencia que son las más altas. Ahora cambiemos la representación gráfica y en lugar que sea una gráfica en tres dimensiones la vamos a crear en 2 dimensiones donde x es la parte real, y es la parte imaginaria y el número de iteraciones es el color del punto

n=900;
nPuntos=n*n;
lienzo = zeros(n, n);
rangox=[-2 2];
rangoy=[-2 2];
xs=linspace(rangox(1),rangox(2),n); %valores reales
ys=linspace(rangoy(1),rangoy(2),n); %valores imaginarios
[x, y] = meshgrid(xs, ys); %combinaciones de cada valor
ks=zeros(n,1);
for j=1:nPuntos
    datos.x0=x(j)+i*y(j); %se crea el valor complejo con i
    [xr,ks(j)]=nr(datos);
    lienzo(j)=ks(j);
end
m = max(lienzo(:));
colores=flipud(bone(m));
colormap(colores); %cambia el mapa de colores
imagesc(xs,ys,lienzo); %crea la imagen en dos dimensiones
axis tight equal

A éste gráfico se le conoce como el Fractal de Newton-Raphson. El gráfico presenta un comportamiento fractal, es decir, si hacemos un acercamiento en cualquier sección se seguirá viendo el mismo comportamiento. Probemos ahora con la siguiente función y su derivada

datos.funcion  = @(x) x.^6+x.^3-1;
datos.derivada = @(x) 6*x.^5+3*x.^2;

n=900;
nPuntos=n*n;
lienzo = zeros(n, n);
rangox=[-2 2];
rangoy=[-2 2];
xs=linspace(rangox(1),rangox(2),n);
ys=linspace(rangoy(1),rangoy(2),n);
[x, y] = meshgrid(xs, ys);
ks=zeros(n,1);
for j=1:nPuntos
    datos.x0=x(j)+i*y(j);
    [xr,ks(j)]=nr(datos);
    lienzo(j)=ks(j);
end
m = max(lienzo(:));
colores=flipud(bone(m));
colormap(colores); %cambia el mapa de colores
imagesc(xs,ys,lienzo); %crea la imagen en dos dimensiones
axis tight equal

Probremos ahora con la siguiente función y su derivada

datos.funcion  = @(x) x.^3-2*x+2;
datos.derivada = @(x) 3*x.^2-2;

n=900;
nPuntos=n*n;
lienzo = zeros(n, n);
rangox=[-2 2];
rangoy=[-2 2];
xs=linspace(rangox(1),rangox(2),n);
ys=linspace(rangoy(1),rangoy(2),n);
[x, y] = meshgrid(xs, ys);
ks=zeros(n,1);
for j=1:nPuntos
    datos.x0=x(j)+i*y(j);
    [xr,ks(j)]=nr(datos);
    lienzo(j)=ks(j);
end
niveles = min(lienzo(:)) + [0:15];
contourf(xs,ys,lienzo,niveles); %crea la imagen en dos dimensiones
colormap(fliplr(colormap('parula')))
colorbar
axis tight equal

Ahora hagamos un acercamiento en el rango [-1,-0.5]

Estos son algunos otros ejemplos de fractales de Newton/Raphson