notes about my code examples here on ComputingSkillSet.com<\/a>.<\/p>\n\n\n\n#!\/usr\/bin\/env python3\n# -*- coding: utf-8 -*-\n\"\"\"\nCreated on Sat Sep 28 05:53:07 2019\n\nNewton-Fractal example code. \nSolves the Newton-method for a given function\non a grid of complex initial guesses.\nProduces a color-coded graph of the fractal.\nIncludes several presets for function, but is easily extendable.\n\n@author: Andreas Krassnigg, ComputingSkillSet.com\/about\/\n\nThis work is licensed under a \nCreative Commons Attribution-ShareAlike 4.0 International License.\nMore information: https:\/\/creativecommons.org\/licenses\/by-sa\/4.0\/\n\n\"\"\"\n\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport datetime\n\n\n# initialize a dictionary of list of the roots\nrootlist = {}\n\n# function definitions to evaluate f\/f' at x. Add your own as needed.\n# each function definition must include a list of the roots of this function\n# root list can be in any order and restricted to finite number\n\n# example 1: polynomial function with four roots\ndef npe1(x):\n return (x**2-1)*(x**2+1)\/(2*x*(x**2-1)+2*x*(x**2+1))\nrootlist['npe1'] = [-1, 1, -1j, 1j]\n\n\n# example 2: function with three roots on the unit circle\ndef npe2(x):\n return (x**3-1)\/(3*x**2)\nrootlist['npe2'] = [-.5-0.8660254037844386j,-.5+0.8660254037844386j, 1]\n\n\n# example 2: function with twelve roots on the unit circle\ndef npe3(x):\n return (x**12-1)\/(12*x**11)\nrootlist['npe3'] = [-.5-0.8660254037844386j,-.5+0.8660254037844386j,.5-0.8660254037844386j,.5+0.8660254037844386j,-.5j-0.8660254037844386,-.5j+0.8660254037844386,.5j-0.8660254037844386,.5j+0.8660254037844386, 1,-1,1.j,-1.j]\n\n\n# example 7: function with four roots, all real\ndef npe7(x):\n return (x+2.)*(x+1.5)*(x-0.5)*(x-2.)\/((x+1.5)*(x-0.5)*(x-2.) + (x+2.)*(x-0.5)*(x-2.) +(x+2.)*(x+1.5)*(x-2.) + (x+2.)*(x+1.5)*(x-0.5)) \nrootlist['npe7'] = [-2, -1.5, 0.5, 2]\n\n\n# example 9: function with four roots, one multiple\ndef npe9(x):\n return (x+2)*(x+1.5)**2*(x-0.5)*(x-2)\/((x+1.5)**2*(x-0.5)*(x-2) + 2*(x+2)*(x+1.5)*(x-0.5)*(x-2) +(x+2)*(x+1.5)**2*(x-2) + (x+2)*(x+1.5)**2*(x-0.5) )\nrootlist['npe9'] = [-2, -1.5, 0.5, 2]\n\n\n# example 10: sine function\ndef npe10(x):\n return np.tan(x)\nrootlist['npe10'] = [0]\nfor i in range(1,10):\n rootlist['npe10'].extend([i*np.pi,-i*np.pi])\n\n\n\n# define a function that can id a root from the rootlist\ndef id_root(zl,rlist):\n findgoal = 1.e-10 * np.ones(len(zl))\n rootid = -1 * np.ones(len(zl))\n for r in rlist:\n # check for closeness to each root in the list\n rootid = np.where(np.abs(zl-r* np.ones(len(zl))) < findgoal, np.ones(len(zl)) * rlist.index(r), rootid)\n \n return rootid\n\n\n\n# define parameters for plotting - adjust these as needed for your function\n \n# define left and right boundaries for plotting \n# for overview plots:\ninterval_left = -2.1\ninterval_right = 2.1\ninterval_down = -2.1\ninterval_up = 2.1\n\n# for detailed plots (adjust as needed):\n#interval_left = 1.15\n#interval_right = 2.\n#interval_down = -0.25\n#interval_up = 0.25\n\n# set number of grid points on x and y axes for plotting \n# use 100 for testing plotting ranges, 1000 for nice plots and 2000 for nicer plots\nnum_x = 1000\nnum_y = 1000\n\n# set desired precision and max number of iterations\n# keep precision goal smaller than findgoal (root matching) above\nprec_goal = 1.e-11\n# max number of iterations. Is being used in a vectorized way. \n# 50 is a good minimal value, sometimes you need 500 or more\nnmax = 200\n\n# timing \nprint('Started computation at '+str(datetime.datetime.now()))\n\n\n# define x and y grids of points for computation and plotting the fractal\nxvals = np.linspace(interval_left, interval_right, num=num_x)\nyvals = np.linspace(interval_down, interval_up, num=num_y)\n\n\n# the following defines function to solve and plot. Jump to end of code to choose function\ndef plot_newton_fractal(func_string, perfom_shading=False):\n \n # create complex list of points from x and y values\n zlist = []\n for x in xvals:\n for y in yvals:\n zlist.append(x + 1j*y)\n \n # initialize the arrays for results, differences, loop counters \n reslist = np.array(zlist)\n reldiff = np.ones(len(reslist))\n counter = np.zeros(len(reslist)).astype(int)\n # initialize overall counter for controlling the while loop\n overallcounter = 0\n # vectorize the precision goal\n prec_goal_list = np.ones(len(reslist)) * prec_goal\n # iterate while precision goal is not met - vectorized\n while np.any(reldiff) > prec_goal and overallcounter < nmax:\n \n # call function as defined above and \n # compute iteration step, new x_i, and relative difference\n diff = eval(func_string+'(reslist)')\n z1list = reslist - diff\n reldiff = np.abs(diff\/reslist)\n \n # reset the iteration\n reslist = z1list\n \n # increase the vectorized counter at each point, or not (if converged)\n counter = counter + np.greater(reldiff, prec_goal_list )\n # increase the control counter\n overallcounter += 1\n \n # get the converged roots matched up with those predefined in the root list\n nroot = id_root(z1list,rootlist[func_string]).astype(int)\n \n # add information about number of iterations to the rood id for shaded plotting\n if perfom_shading == True:\n nroot = nroot - 0.99*np.log(counter\/np.max(counter))\n\n # uncomment those in case of doubt\n# print(reslist)\n# print(counter)\n# print(nroot)\n \n # get the data into the proper shape for plotting with matplotlib.pyplot.matshow\n nroot_contour = np.transpose(np.reshape(nroot,(num_x,num_y)))\n \n # timing to see difference in time used between calculation and plotting\n print('Finished computation at '+str(datetime.datetime.now()))\n \n # create an imshow plot \n plt.figure()\n \n #label the axes\n plt.xlabel(\"$Re(z)$\", fontsize=16)\n plt.ylabel(\"$Im(z)$\", fontsize=16)\n \n # plots the matrix of data in the current figure. Interpolation isn't wanted here.\n # Change color map (cmap) for various nice looks of your fractal\n plt.matshow(nroot_contour, fignum=0, interpolation='none', origin='lower', cmap='hot')\n \n # remove ticks and tick labels from the figure\n plt.xticks([])\n plt.yticks([])\n \n # save a file of you plot. 200 dpi is good to match 1000 plot points. \n # Increasing dpi produces more pixels, but not more plotting points.\n # Number of plotting points and pixels are best matched by inspection.\n plt.savefig('newton-fractal-example-plot-'+func_string+'.jpg', bbox_inches='tight', dpi=200)\n \n plt.close()\n # timing step\n print('Finished creating matshow plot at '+str(datetime.datetime.now()))\n\n\n# call the solution function for func of your choice\n# also, switch shading via number of iterations on or off\nplot_newton_fractal('npe1', perfom_shading=True)\n\n# final timing check\nprint('Finished computation and plotting at '+str(datetime.datetime.now()))\n\n<\/pre>\n\n\n\n