by Anand | tags: [ fractal python newton newton-raphson calculus mathematics visualization matplotlib numpy ]
Fractals with Python - The Newton Set
This is the fourth post of the fractal series. Today I am going to iterate to a new fractal family, namely the Newton fractal set.
The Mathematics
The Newton-Raphson technique is a well-known iterative technique for finding roots of a function by continually improving an approximation to the root using the function derivative.
If you have a function $f(x)$ and an approximation to its root $a$, then this technique calculates the next better approximation to the root as,
$$a^′ = a - f(a)/f^{′}(a)$$
This step is repeated till successive values converge quickly to the root.
In effect, you are drawing a tangent to the curve at the point of the original approximation and following the tangent to where it intersects the x-axis.
If one starts with a bad approximation however, one may end up with roots that diverge instead of converge as in the case where one picks a point at the top of a parabola that cross the x-axis at two places. It may end up spanning the entire real line while oscillating back and forth.
This kind of behaviour where a small area expands to a larger area and back and forth leads to self-similar shapes - in other words, Fractals.
Hence, if we start a Newton-Raphson iteration at each point on the complex plane, run it until it converges to a toleration close to a root and then color the starting point according to either the root or the number of iterations, we get a set of fractal shapes.
Let us see how this can be implemented in Python code.
Newton Set for $f(z) = z^3 - 1$
For this function, one can compute the (complex) roots as,
- $z =1$
- $z = -0.5 + \frac{\sqrt{3}}{2}j$
- $z = -0.5 - \frac{\sqrt{3}}{2}j$
We will create the code step by step. First lets define the Python function. It will be defined specifically for the given polynomial function $f(z) = z^3 - 1$.
def newton_set(width, height, zoom=1, x_off=0, y_off=0, niter=256):
""" Fractals using newton-raphson """
# Pixels array
pixels = np.arange(width*height*3, dtype=np.uint32).reshape(height, width, 3)
Let’s define the drawing area for this in the square defined by the vertices $$(-2.5, -2.5), (-2.5, 2.5), (2.5, 2.5), (2.5, -2.5)$$.
# drawing area
xa, xb, ya, yb = -2.5, 2.5, -2.5, 2.5
The step size for the derivative, the epsilon (max error) are defined next.
h = 1e-7 # step size for numerical derivative
eps = 1e-3 # max error allowed
Since we are going to color using roots, let’s define the roots for this function, along with a color multiplication factor.
# Bounding roots
r1 = 1
r2 = complex(-0.5, math.sin(2*math.pi/3))
r3 = complex(-0.5, -1*math.sin(2*math.pi/3))
# Color multiplication factor
# NOTE: Increasing this darkens the image and lightens it otherwise
multcol = 5
Once the definitions are out of the way its time to implement the main iterative loop of the function having the fractal computation logic.
for y in range(height):
zy = (y + y_off) * (yb - ya) / (zoom*(height - 1)) + ya
for x in range(width):
zx = (x + x_off) * (xb - xa) / (zoom*(width - 1)) + xa
# Mapping real -> complex plane
z = complex(zx, zy)
count = 0
for i in range(niter):
# complex numerical derivative
dz = (fcube(z + complex(h, h)) - fcube(z)) / complex(h, h)
if dz == 0:
break
count += 1
if count > 255:
break
znext = z - fcube(z) / dz # Newton iteration
if abs(znext - z) < eps: # stop when close enough to any root
break
z = znext
# Pixels colored using the roots
if abs(z-r1)<eps:
# color red
pixels[y,x] = (255 - count*multcol, 0, 0)
elif abs(z-r2)<=eps:
# color green
pixels[y,x] = (0, 255 - count*multcol, 0)
elif abs(z-r3)<=eps:
# color blue
pixels[y,x] = (0, 0, 255 - count*multcol)
Last but not the least, is the definition for $fcube$ .
def fcube(z):
return z ** 3 - 1
The code is well commented and is hopefully quite self-explanatory.
Briefly, we take a point (pixel) in the screen, map it to the complex plane as $z$ and calculate $f′(z)$ at the point. The next point is calculated using the newton-rapshon formula. When the point reaches close to the root, using an error value (epsilon), we stop the iteration.
The color value for the pixel is calculated using the roots by determining how close the point is to a given root.
The function will be called as usual by the display
wrapper function.
def display(width=1024, height=1024, niter=1024, zoom=1, x_off=0, y_off=0):
""" Display a newton-raphson fractal """
pimg = newton_set(width, height, zoom=zoom, x_off=x_off, y_off=y_off, niter=niter)
plt.axis('off')
plt.imshow(pimg)
plt.show()
Visualizing the Set
The function defined earlier produces the following fractal image.
Fractal image for $f(z) = z^3 - 1$ colored using roots
Let’s zoom into this a little bit. Using,
>>> newton.display(zoom=3, x_off=-50, y_off=1000)
It produces an image like the one below with the fractal zoomed in, showing the self-similarity and repeating patterns.
Fractal image for $f(z) = z^3 - 1$ colored using roots (at zoom = 3, y_offset=1000, x_offset=-50)
Generalizing the Fractal Code
The issue with the fractal code as implemented is that for each function we need to find its roots and create custom coloring code for coloring each pixel according to the boundaries defined by the roots. A more general function can be written by coloring the fracal according to the number of iterations than by using the roots.
While doing this, we will also conveniently parallelize our code to take advantage of the CPU cores using the techniques described in previous articles.
Here is the code that calculates each row of the fractal.
def newton_set_calc_row(y, width, height, function, niter=256, x_off=0, y_off=0, zoom=1):
""" Calculate one row of the newton set with size width x height """
row_pixels = np.arange(width, dtype=np.uint32)
# drawing area
xa, xb, ya, yb = -2.5, 2.5, -2.5, 2.5
zy = (y + y_off) * (yb - ya) / (zoom*(height - 1)) + ya
h = 1e-7 # step size for numerical derivative
eps = 1e-3 # max error allowed
for x in range(width):
# calculate the initial real and imaginary part of z,
# based on the pixel location and zoom and position values
zx = (x + x_off) * (xb - xa) / (zoom*(width - 1)) + xa
z = complex(zx, zy)
for i in range(niter):
# complex numerical derivative
dz = (function(z + complex(h, h)) - function(z)) / complex(h, h)
if dz == 0:
break
znext = z - a*function(z) / dz # Newton iteration
if abs(znext - z) < eps: # stop when close enough to any root
break
z = znext
# Color according to iteration count
rgb = (i % 16 * 32, i % 8 * 64, i % 4 * 64)
row_pixels[x] = rgb2int(rgb)
return y,row_pixels
Here are the functions that converts a $(r,g,b)$ tuple to an int and viceverza.
def rgb2int(rgb):
""" Convert (r,g,b) tuple into an int """
red = rgb[0]
green = rgb[1]
blue = rgb[2]
return (red << 16) + (green << 8) + blue
def int2rgb(rgbint):
blue = rgbint & 255
green = (rgbint >> 8) & 255
red = (rgbint >> 16) & 255
return (red, green, blue)
And here is the calling function which performs parallel computation across all CPU cores by using the above function.
import functools
import multiprocessing as mp
def newton_set_mp(width, height, function, zoom=1, x_off=0, y_off=0, niter=256):
""" Newton-raphson fractal set with multiprocessing """
w,h = width, height
pixels = np.arange(w*h*3, dtype=np.uint32).reshape(h, w, 3)
# print('Starting calculation using',width, height,cx,cy)
pool = mp.Pool(mp.cpu_count())
newton_partial = functools.partial(newton_set_calc_row,
width=width,height=height, function=function,
niter=niter,zoom=zoom,x_off=x_off,y_off=y_off)
for y,row_pixel in pool.map(newton_partial, range(h)):
for x in range(w):
pixels[y, x] = np.array(int2rgb(row_pixel[x]))
return pixels
Finally, here is the new display
function.
def display(function, width=1024, height=1024, niter=1024, zoom=1, x_off=0, y_off=0):
""" Display a newton-raphson fractal """
pimg = newton_set_mp(width, height, function, zoom=zoom,x_off=x_off, y_off=y_off, niter=niter)
plt.axis('off')
plt.imshow(pimg)
plt.show()
NOTE: Observe how the `display` function accepts an additional parameter `function`. This parametrizes it to accept any function and plot its fractal.
I will not discuss the implementation in detail since the coding techniques used are same as those used for the Julia and Mandelbrot fractals in previous iterations of this series.
Here is the fractal where for function = fcube
.
Fractal image for $f(z) = z^3 - 1$ colored using iterations
Fractal Images for Polynomial Functions
Since our code is parametrized to accept any function, we can experiment with it to generate images for multiple functions.
Let us define a few functions in higher degrees of z
and plot their fractals.
def fquad(z):
return z**4 - 1
def fsix(z):
return z**6 + z**3 - 1
def f8(z):
return z**8 + 15*z**4 - 16
Here are their respective fractal images.
Fractal image for $f(z) = z^4 - 1$
Fractal image for $f(z) = z^6 + z^3 - 1$
Fractal image for $f(z) = z^8 + 15z^4 - 16$
There are some interesting features about these images. For example, you may have noted that the number of $arms$ for each image exactly matches the degree of its function.
Fractal Images - Trigonometric Functions
Here are some images the code generates for some well known trigonometric functions.
def sin(z):
return math.sin(abs(z))
Fractal image for $f(z) = sin(z)$
def cos(z):
return math.cos(abs(z))
Fractal image for $f(z) = cos(z)$
def tan(z):
return math.tan(abs(z))
Fractal image for $f(z) = tan(z)$
Experiments and Follow-Up
One can experiment with the color scheme used in the code for coloring using iterations and getting various color $palettes$ . It is also suggested as an exercise to the reader to modify the first function - that uses the roots to color the image - and extend it for a higher degree function say $f(x) = x^4 - 1$.
For those with a more mathematical inclination, it may be interesting to try and explain the $shape$ of these fractal images in relation to the generating function and also try and play around with more complex functions.
If you liked the article and want to inspect the code in depth, you can heade over to the code for this in our repository.
This was hopefully as entertaining and educational to the reader as it was to me in writing it. I will come back in the future with another edition of this series. Till then enjoy your coffee, your code and your fractals!
References
Note that name and e-mail are required for posting comments