Consider the function $$ f(x) = x^5 + 1 $$ The roots lie on the unit circle in the complex plane. We can write $$ -1 = \exp(i(2n+1)\pi), \qquad n=0,1,2,\ldots $$ The fifth roots are $$ (-1)^{1/5} = \exp(i(2n+1)\pi/5), \qquad n=0,1,2,3,4 $$
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
from numpy import abs,exp,array,real,imag,pi
import matplotlib.pyplot as plt
Lets compute the exact roots and plot them on a graph.
n = array([0,1,2,3,4])
r = exp(1j*(2*n+1)*pi/5) # 5'th roots of -1
plt.plot(real(r),imag(r),'o')
plt.grid(True)
plt.axis('equal')
for x in r:
print(x,x**5)
(0.8090169943749475+0.5877852522924731j) (-1.0000000000000002+0j) (-0.30901699437494734+0.9510565162951536j) (-1+4.440892098500626e-16j) (-1+1.2246467991473532e-16j) (-1+6.123233995736766e-16j) (-0.30901699437494756-0.9510565162951535j) (-1+7.216449660063518e-16j) (0.8090169943749473-0.5877852522924732j) (-1+6.661338147750939e-16j)
Define the function and its derivative.
def f(x):
return x**5 + 1.0
def df(x):
return 5.0*x**4
Now implement Newton method as a function.
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
We call Newton method with a complex initial guess for the root.
x0 = 1.0 + 1.0j
x = newton(f,df,x0)
print('Root = ',x)
0 (0.85+0.8j) 1.4844926068747277 1 (0.7869450486819533+0.6530228729743794j) 0.3587868496977192 2 (0.8000098256925355+0.5887168663596131j) 0.04466991909204323 3 (0.8091271135998733+0.587660739881959j) 0.0008311338875735857 4 (0.8090169566882729+0.5877852118808937j) 2.762870121656608e-07 5 (0.8090169943749507+0.5877852522924782j) 3.014107274014632e-14 6 (0.8090169943749475+0.5877852522924731j) 2.220446049250313e-16 7 (0.8090169943749475+0.5877852522924731j) 2.220446049250313e-16 Root = (0.8090169943749475+0.5877852522924731j)