The Julia set

Question P7.2.4

The Julia set associated with the complex function $f(z) = z^2 + c$ may be depicted using the following algorithm.$\renewcommand\Re{\operatorname{Re}}\renewcommand\Im{\operatorname{Im}}$

For each point, $z_0$, in the complex plane such that $-1.5 \le \Re[z_0] \le 1.5$ and $-1.5 \le \Im[z_0] \le 1.5$, iterate according to $z_{n+1} = z_n^2 + c$ where $c$ is a some (complex) constant. Colour the pixel in an image corresponding to this region of the complex plane according to the number of iterations required for $|z|$ to exceed some critical value, $|z|_\mathrm{max}$ (or black if this does not happen before a certain maxmimum number of iterations $n_\mathrm{max}$).

Write a program to plot the Julia set for $c = -0.1 + 0.65j$, using $|z|_\mathrm{max}=10$ and $n_\mathrm{max}=500$.

Solution

Click here for a solution

The code below is one way of plotting the required Julia set. There are many different ways of picking the pixel colours.

importnumpyasnpimportmatplotlib.pyplotaspltimportmatplotlib.cmascm# Image width and height; parameters for the plotim_width,im_height=500,500c=complex(-0.1,0.65)zabs_max=10nit_max=1000xmin,xmax=-1.5,1.5xwidth=xmax-xminymin,ymax=-1.5,1.5yheight=ymax-yminjulia=np.zeros((im_width,im_height))forixinrange(im_width):foriyinrange(im_height):nit=0# Map pixel position to a point in the complex planez=complex(ix/im_width*xwidth+xmin,iy/im_height*yheight+ymin)# Do the iterationswhileabs(z)<=zabs_maxandnit<nit_max:z=z**2+cnit+=1shade=1-np.sqrt(nit/nit_max)ratio=nit/nit_maxjulia[ix,iy]=ratiofig,ax=plt.subplots()ax.imshow(julia,interpolation='nearest',cmap=cm.hot)# Set the tick labels to the coordinates of z0 in the complex planextick_labels=np.linspace(xmin,xmax,xwidth/0.5)ax.set_xticks([(x-xmin)/xwidth*im_widthforxinxtick_labels])ax.set_xticklabels(['{:.1f}'.format(xtick)forxtickinxtick_labels])ytick_labels=np.linspace(ymin,ymax,yheight/0.5)ax.set_yticks([(y-ymin)/yheight*im_heightforyinytick_labels])ax.set_yticklabels(['{:.1f}'.format(ytick)forytickinytick_labels])plt.show()