Consider the function $$ f(x) = x^5 + 1 $$ The roots lie on the unit circle in the complex plane.
from numpy import abs
def f(x):
return x**5 + 1.0
def df(x):
return 5.0*x**4
def newton(f,df,x0,M=100,eps=1.0e-15):
x = x0
for i in range(M):
dx = - f(x)/df(x)
x = x + dx
print i,x,abs(f(x))
if abs(dx) < eps * abs(x):
return x
x0 = 1.0 + 1.0j
x = newton(f,df,x0)
0 (0.85+0.8j) 1.48449260687 1 (0.786945048682+0.653022872974j) 0.358786849698 2 (0.800009825693+0.58871686636j) 0.044669919092 3 (0.8091271136+0.587660739882j) 0.000831133887574 4 (0.809016956688+0.587785211881j) 2.76287012166e-07 5 (0.809016994375+0.587785252292j) 3.01410727401e-14 6 (0.809016994375+0.587785252292j) 2.22044604925e-16 7 (0.809016994375+0.587785252292j) 2.22044604925e-16