In [2]:
from numpy import *
import matplotlib.pyplot as plt
%matplotlib inline
In [10]:
nr = 5
nc = 6
r = 1.5

x = linspace(-r,r,nc)
y = linspace(-r,r,nr)

x
Out[10]:
array([-1.5, -0.9, -0.3,  0.3,  0.9,  1.5])
In [11]:
y
Out[11]:
array([-1.5 , -0.75,  0.  ,  0.75,  1.5 ])
In [12]:
X,Y = meshgrid(x,y)
In [13]:
Z = X + 1j*Y
In [14]:
Z
Out[14]:
array([[-1.5-1.5j , -0.9-1.5j , -0.3-1.5j ,  0.3-1.5j ,  0.9-1.5j ,
         1.5-1.5j ],
       [-1.5-0.75j, -0.9-0.75j, -0.3-0.75j,  0.3-0.75j,  0.9-0.75j,
         1.5-0.75j],
       [-1.5+0.j  , -0.9+0.j  , -0.3+0.j  ,  0.3+0.j  ,  0.9+0.j  ,
         1.5+0.j  ],
       [-1.5+0.75j, -0.9+0.75j, -0.3+0.75j,  0.3+0.75j,  0.9+0.75j,
         1.5+0.75j],
       [-1.5+1.5j , -0.9+1.5j , -0.3+1.5j ,  0.3+1.5j ,  0.9+1.5j ,
         1.5+1.5j ]])
In [15]:
def g(z): return 2/3*z + 1/3/z**2
In [16]:
g(Z)
Out[16]:
array([[-1.00000000-1.07407407j, -0.65126233-1.09611688j,
        -0.33149244-1.05478852j,  0.06850756-0.94521148j,
         0.54873767-0.90388312j,  1.00000000-0.92592593j],
       [-0.92888889-0.59481481j, -0.55620452-0.73888441j,
        -0.56992998-0.85231426j, -0.16992998-0.14768574j,
         0.64379548-0.26111559j,  1.07111111-0.40518519j],
       [-0.85185185+0.j        , -0.18847737+0.j        ,
         3.50370370+0.j        ,  3.90370370+0.j        ,
         1.01152263+0.j        ,  1.14814815+0.j        ],
       [-0.92888889+0.59481481j, -0.55620452+0.73888441j,
        -0.56992998+0.85231426j, -0.16992998+0.14768574j,
         0.64379548+0.26111559j,  1.07111111+0.40518519j],
       [-1.00000000+1.07407407j, -0.65126233+1.09611688j,
        -0.33149244+1.05478852j,  0.06850756+0.94521148j,
         0.54873767+0.90388312j,  1.00000000+0.92592593j]])
In [29]:
nr = 1000
nc = 1000
r = 1.5

x = linspace(-r, r,nc)
y = linspace( r,-r,nr)
X,Y = meshgrid(x,y)

Z = X + 1j*Y

def g(z): return 2/3*z + 1/3/z**2   # Newton map for f(z) = z^3 - 1

# create a color image
a = zeros((nr,nc,3),dtype=float)#uint8)
# for convenience, set up views into each color layer
red   = a[:,:,0]
green = a[:,:,1]
blue  = a[:,:,2]

redroot   = 1.
greenroot = -1/2 + 1j*sqrt(3)/2
blueroot  = -1/2 - 1j*sqrt(3)/2

for k in range(20):
    # add 1 to red layer at pixels where Z now close to red root, etc.
    # use boolean indexing
    red  [ abs(Z-  redroot) < .1    ] += 1    # add one to every pixel where Z now close to red root
    green[ abs(Z-greenroot) < .1    ] += 1    # add one to every pixel where Z now close to green root
    blue [ abs(Z- blueroot) < .1    ] += 1    # add one to every pixel where Z now close to blue root

    Z = g(Z)

from scipy.misc import imsave, imshow
imsave('newtonconvergence.png',a)
print(a.max())
plt.figure(figsize=(6,6))
plt.imshow(a,interpolation='none',extent=[-r,r,-r,r]);
20.0
In [ ]:
Why is my imshow not showing gradations of color??