Diving Into the Deluge of Data :: Lab 3 :: Numerical Methods: Roots and Fractals

Lab 3: Numerical Methods: Roots and Fractals

This labs explores algorithmic techniques for successive approximations to the square root function. This includes the bisection method and Newton's method. Moving beyond square roots to the complex cube roots of unity, this lab also explores ways of visualizing common areas of convergence resulting in a well-known Julia set called the Newton basin.

Step 0: Lab Preparation

Step 1: Source Code

Step 2: Calculating Square Roots using Successive Approximations

The file sqrt.py contains function headers for sqrt_newton and sqrt_bisect. Implement these functions using x/2 as the starting value in sqrt_newton and the values 0 and x as the starting low and high values respectively in sqrt_bisect. 

    def sqrt_newton(x, iter):
        """
        Use Newton's Method for 'iter' iterations
        to approximate sqrt(x)
        """
    def sqrt_bisect(x, iter):
        """
        Use bisection / binary search method for 'iter' iterations
        to approximate sqrt(x)
        """

You can always test your functions out by running python in interpreter mode and typing 

	>>> import sqrt
	>>> sqrt.sqrt_newton(25,3)
	5.011394106532552
	>>> sqrt.sqrt_bisect(25,3)
	4.6875

Implement the functions absolute_error, list_of_errors, and save_plot. You will find the matplotlib function plt.savefig("somefile.png") useful for saving the current plot to disk. Below you will also find a picture of the plot that you can use as a guide when coding the save_plot function. 

    def absolute_error(square, iterations, sqrt_fn):
      """
      Calculate the absolute error between the value reported by
      'sqrt_fn' and the actual square root of 'square'.  Absolute error is the absolute
      difference between the correct square root and the value given by 'sqrt_fn' after
      'iterations'"""
    def list_of_errors(square, max_iterations, sqrt_fn):
      """
      Return a list of 'max_iterations' absolute error values for 'sqrt_fn' where 
      the error value at index i represents the absolute error between 'sqrt_fn' 
      on 'square' after i iterations and the actual square root.
      """
    def save_plot(newton_err, bisect_err, filename):
      """
      Save a figure to 'filename' that shows absolute error with respect to 
      number of iterations for both Newton's method and the bisection method.
      """

Calling 

$ python sqrt.py 25 15 error.png

from the command line should produce a file called error.png that you can open in Preview by typing $ open error.png. It should look like this:

Step 3: Making Images

Included in this week's source code is a module called image.py that provides a simple interface for creating images, drawing points, and saving images. 
    from PIL import Image, ImageDraw

    def create_image(width, height):
      return Image.new("RGB", (width, height), (255, 255, 255))

    def draw_point(image, x, y, color):
      ImageDraw.Draw(image).point((x,y), color)

    def save_image(image, filename):
      image.save(filename, "PNG")

We can use the above interface to create images. For example, here's a 100x100 pixel image with a 20-pixel wide horizontal stripe. Note that the coordinate system for images runs from top-left (0,0) to bottom-right (100,100). Colors are just 3-tuples of RGB (Red, Green, and Blue) values that each range from 0 to 255. 

      >>> im = image.create_image(100,100)
      >>> for x in range(100):
      ...   for y in range(40,60):
      ...     image.draw_point(im, x, y, (255, 0, 0))
      ... 
      >>> image.save_image(im, "redstripe.png")

Step 4: Newton's Basin

If z is a complex number then the function f(z)=z3-1 has three different complex roots at

To which root Newton's method converges depends entirely on the point at which it starts. Surprisingly, this starting point is very sensitive. In fact, if one visualizes the complex plane so that one colors each complex point according to the root on which it converges using Newton's method, you obtain the following fractal known as Newton's fractal

Recall that for some polynomial f(z), Newton's method finds successive approximations of a root of f(z) by letting the approximation zn at iteration n be given by the following assignment:

zn = zn-1 - f(zn-1) / f'(zn-1)

so for f(z)=z3-1 this translates to

zn = zn-1 - (zn-13-1) / 3(zn-12)

Examine the source code in basin.py. One function closest_root, is defined for you. Three others you must define.

You should interpret the interpolate_value functions as taking a value x in the range [0,length-1] and returning the corresponding position in the range (y1,y2). That is you should return y1 + (x / length * (y2 - y1))

In particular, the draw function should do the following: 

    def cube_root_of_unity(z, iter):
        """Use Newtown's Method for 'iter' iterations
           to find a cube root of unity.  Start the approximation
           at point z"""
    def interpolate_value(x,length,y1,y2):
        """Linearly interpolate the value x in the range [0,length-1] 
           to range [y1,y2]"""
    def draw(output):
        """
        1.  Create an image with dimesions IMAGE_WIDTH x IMAGE_HEIGHT
        2.  For each pixel (i,j) in the image, linearly interpolate 
            that pixel into the complex plane.  The top-left corner
            of the image should correspond to the bttom-left corner 
            of the complex plane.
        3.  If the pixel gets interpolated to (0,0) then skip it.
        4.  Given the complex point, check to which root it converges 
            using Newton's method.  Color the pixel
            (a) RED if it's ZROOT1,
            (b) BLUE if it's ZROOT2, and
            (c) GREEN if it's ZROOT3
        5.  Save the resuling image as 'filename'"""

You should be able to run 

$ python basin.py basin.png
and produce the image below.