Newton's Method
system:sage

{{{id=0|
f(x)=x^3+2*x^2+3*x+4
///
}}}

{{{id=1|
z={}
z[0]=3+7*I
z[1]=z[0]-f(z[0])/f.derivative(x)(z[0])
print z[1]
print CDF(z[1])
z[2]=z[1]-f(z[1])/f.derivative(x)(z[1])
print z[2]
print CDF(z[2])
///


23424/4963*I + 8752/4963
1.7634495265 + 4.71972597219*I
2218930600851456/686077676398303*I + 643334213751084/686077676398303
0.937698799833 + 3.23422649823*I
}}}

{{{id=3|
for i in range(10):
    z[i+1]=z[i]-f(z[i])/f.derivative(x)(z[i])
    print CDF(z[i+1])
///


1.7634495265 + 4.71972597219*I
0.937698799833 + 3.23422649823*I
0.392132013081 + 2.30429653276*I
0.0500091493907 + 1.78561924162*I
-0.125366996437 + 1.57918320368*I
-0.172268166359 + 1.54704061589*I
-0.174682973653 + 1.54686531921*I
-0.174685404294 + 1.54686888723*I
-0.17468540428 + 1.54686888723*I
-0.17468540428 + 1.54686888723*I
}}}

{{{id=7|
def newton(func, guess=0, acc=.00001):
    curr_guess = CDF(guess)
    delta = CDF(1+acc)
    count = 0
    deriv = func.derivative(x)
    while delta.abs() > acc:
        count += 1
        if deriv(curr_guess)==0:
            return curr_guess #Not really the right thing to do.
        delta = func(curr_guess)/deriv(curr_guess)
        curr_guess = curr_guess - delta
        if count>10000:
            raise Exception("Took too long.")
    return curr_guess
///
}}}

{{{id=8|
newton(g, guess=0)
///

0
}}}

{{{id=9|
g(x)=x^3-1
///
}}}

{{{id=10|
newton(g, 0)
///


Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/lipshitz/.sage/sage_notebook/worksheets/admin/2/code/110.py", line 7, in <module>
    exec compile(ur'newton(g, _sage_const_0 )' + '\n', '', 'single')
  File "", line 1, in <module>
    
  File "/Users/lipshitz/.sage/sage_notebook/worksheets/admin/2/code/107.py", line 14, in newton
    delta = func(curr_guess)/deriv(curr_guess)
  File "element.pyx", line 1271, in sage.structure.element.RingElement.__div__ (sage/structure/element.c:10461)
  File "expression.pyx", line 1863, in sage.symbolic.expression.Expression._div_ (sage/symbolic/expression.cpp:10771)
ZeroDivisionError: Symbolic division by zero
}}}

{{{id=11|
def which_sol(guess):
    "Figures out which solution Newton's method converges to for the equation z^3=1"
    g(x)=x^3-1
    sol=newton(g(x),guess)
    if CDF(sol-1).abs() < .1: return 1
    if CDF(sol+1/2-sqrt(3)*I/2).abs()<.1: return 2
    if CDF(sol+1/2+sqrt(3)*I/2).abs()<.1: return 3
    return 4
///
}}}

{{{id=12|
print which_sol(-2+2*I)
///


2
}}}

{{{id=16|
size=100
M=matrix(size,size)
for i in range(size):
    for j in range(size):
        M[i,j]=which_sol(i*2/size+j*2/size*I-1.0001-I)
///
}}}

{{{id=19|
matrix_plot(M,cmap='autumn_r')
///


<html><font color='black'><img src='cell://sage0.png'></font></html>
}}}

{{{id=20|
import matplotlib.cm
matplotlib.cm.datad.keys()
///

['Spectral', 'summer', 'RdBu', 'gist_earth', 'Set1', 'Set2', 'Set3', 'Dark2', 'hot', 'PuOr_r', 'PuBuGn_r', 'RdPu', 'gist_ncar_r', 'gist_yarg_r', 'Dark2_r', 'YlGnBu', 'RdYlBu', 'hot_r', 'gist_rainbow_r', 'gist_stern', 'cool_r', 'cool', 'gray', 'copper_r', 'Greens_r', 'GnBu', 'gist_ncar', 'spring_r', 'gist_rainbow', 'RdYlBu_r', 'gist_heat_r', 'OrRd_r', 'bone', 'gist_stern_r', 'RdYlGn', 'Pastel2_r', 'spring', 'Accent', 'YlOrRd_r', 'Set2_r', 'PuBu', 'RdGy_r', 'spectral', 'flag_r', 'jet_r', 'RdPu_r', 'gist_yarg', 'BuGn', 'Paired_r', 'hsv_r', 'YlOrRd', 'Greens', 'PRGn', 'gist_heat', 'spectral_r', 'Paired', 'hsv', 'Oranges_r', 'prism_r', 'Pastel2', 'Pastel1_r', 'Pastel1', 'gray_r', 'PuRd_r', 'Spectral_r', 'BuGn_r', 'YlGnBu_r', 'copper', 'gist_earth_r', 'Set3_r', 'OrRd', 'PuBu_r', 'winter_r', 'jet', 'bone_r', 'BuPu', 'Oranges', 'RdYlGn_r', 'PiYG', 'YlGn', 'binary_r', 'gist_gray_r', 'BuPu_r', 'gist_gray', 'flag', 'RdBu_r', 'BrBG', 'Reds', 'summer_r', 'GnBu_r', 'BrBG_r', 'Reds_r', 'RdGy', 'PuRd', 'Accent_r', 'Blues', 'Greys', 'autumn', 'PRGn_r', 'Greys_r', 'pink', 'binary', 'winter', 'pink_r', 'prism', 'YlOrBr', 'Purples_r', 'PiYG_r', 'YlGn_r', 'Blues_r', 'YlOrBr_r', 'Purples', 'autumn_r', 'Set1_r', 'PuOr', 'PuBuGn']
}}}

{{{id=21|
size=100
N=matrix(RR,size,size)
for i in range(size):
    for j in range(size):
        new_sol=newton(g,guess=i*2/size+j*2/size*I-1-I)
        N[i,j]=new_sol.real()+new_sol.imag()
///
}}}

{{{id=22|
matrix_plot(N,cmap='autumn_r')
///

<html><font color='black'><img src='cell://sage0.png'></font></html>
}}}

{{{id=23|
size=150
N=matrix(RR,size,size)
for i in range(size):
    for j in range(size):
        new_sol=newton(f,guess=i*4/size+j*4/size*I-2-2*I)
        N[i,j]=new_sol.real()+new_sol.imag()
///
}}}

{{{id=25|
matrix_plot(N,cmap='autumn_r')
///

<html><font color='black'><img src='cell://sage0.png'></font></html>
}}}

{{{id=26|

///
}}}