Welcome, guest | Sign In | My Account | Store | Cart

The code generates complex polynomials that has random real coefficients; each +1 or -1. Later it plots the roots.

Warning: The calculation may take 15 minutes or so!

Python, 66 lines
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66``` ```# Complex Polynomial Roots Fractal # See: Scientific American Magazine, May 2003, "A Digital Slice of Pi" # FB - 201109051 import math import random import time from PIL import Image imgx = 512 imgy = 512 image = Image.new("RGB", (imgx, imgy)) # drawing area xa = -2.0 xb = 2.0 ya = -2.0 yb = 2.0 nmax = 18 # max polynomial degree allowed pmax = 20000 # max number of polynomials to find roots for maxIt = 1000 # cutoff for iteration def polynomial(z, coefficients): # Horner scheme t = complex(0, 0) for c in reversed(coefficients): t = t * z + c return t # Durand-Kerner method for solving polynomial equations eps = 1e-7 # max error allowed def DurandKerner(coefficients): n = len(coefficients) - 1 roots = [complex(0, 0)] * n for k in range(n): roots[k] = complex(random.random(), random.random()) rootsNew = roots[:] flag = True it = 0 # number of iterations while flag: flag = False for k in range(n): temp = complex(1, 0) for j in range(n): if j != k: temp *= roots[k] - roots[j] rootsNew[k] = roots[k] - polynomial(roots[k], coefficients) / temp if abs(roots[k] - rootsNew[k]) > eps: flag = True roots = rootsNew[:] it += 1 if it > maxIt: flag = False return (roots, it) st = time.time() for p in range(pmax): pd = random.randint(2, nmax) coefficients = [complex(0, 0)] * (pd + 1) for k in range(pd): coefficients[k] = complex(math.copysign(1.0, random.randint(0, 1) * 2 - 1), 0) (roots, i) = DurandKerner(coefficients) for k in range(len(roots)): kx = (imgx - 1) * (roots[k].real - xa) / (xb - xa) ky = (imgy - 1) * (roots[k].imag - ya) / (yb - ya) if kx >= 0 and kx <= imgx -1 and ky >= 0 and ky <= imgy - 1: image.putpixel((int(kx), int(ky)), (i % 8 * 32, i % 4 * 64, i % 16 * 16)) image.save("ComplexPolynomialRootsFractal.png", "PNG") print "Duration in seconds: " + str(int(time.time() - st)) ``` FB36 (author) 10 years, 8 months ago

This version uses NumPy library to calculate polynomial roots and it is a lot faster:

``````# Complex Polynomial Roots Fractal
# See: Scientific American Magazine, May 2003, "A Digital Slice of Pi"
# FB - 201109106
import numpy
import math
import random
import time
from PIL import Image
imgx = 512
imgy = 512
image = Image.new("RGB", (imgx, imgy))
# drawing area
xa = -1.5
xb = 1.5
ya = -1.5
yb = 1.5
pdmax = 18 # max polynomial degree allowed
pmax = 50000 # max number of polynomials to find roots for
st = time.time()
for p in range(pmax):
pd = random.randint(2, pdmax) # polynomial degree
coefficients = [random.choice([1, -1]) for k in range(pd + 1)]
roots = numpy.roots(coefficients)
for k in range(len(roots)):
kx = (imgx - 1) * (roots[k].real - xa) / (xb - xa)
ky = (imgy - 1) * (roots[k].imag - ya) / (yb - ya)
if kx >= 0 and kx <= imgx -1 and ky >= 0 and ky <= imgy - 1:
image.putpixel((int(kx), int(ky)), (0, 255, 0))
image.save("ComplexPolynomialRootsFractal.png", "PNG")
print "Calculation time in seconds: " + str(int(time.time() - st))
`````` FB36 (author) 10 years, 8 months ago

Unlike Mandelbrot fractal for example, if you want to draw (zoom) a particular sub-section of this fractal, still the whole fractal must be calculated. Because there is no way to predict which sets of polynomial coefficients would produce roots in the targeted area.

Unless there is a solution for this problem: Given one of its roots (a complex number) is there any polynomial w/ degree less than N that has only a set of coefficients each +1 or -1? (If there is a solution then you could just scan all points in the target area and test each point to decide which ones needs to be plotted.) Created by FB36 on Mon, 5 Sep 2011 (MIT)