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,

  1. $z =1$
  2. $z = -0.5 + \frac{\sqrt{3}}{2}j$
  3. $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.

Newton Fractal Set

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.

Newton Fractal Set

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.

Newton Fractal Set colored using iterations

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.

Newton Fractal Set

Fractal image for $f(z) = z^4 - 1$

Newton Fractal Set

Fractal image for $f(z) = z^6 + z^3 - 1$

Newton Fractal Set

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))
Newton Fractal Set

Fractal image for $f(z) = sin(z)$

def cos(z):
    return math.cos(abs(z))
Newton Fractal Set

Fractal image for $f(z) = cos(z)$

def tan(z):
    return math.tan(abs(z))
Newton Fractal Set

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

  1. https://www.chiark.greenend.org.uk/~sgtatham/newton/

Note that name and e-mail are required for posting comments