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))

2 comments

FB36 (author) 10 years, 3 months ago  # | flag

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, 2 months ago  # | flag

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)
Python recipes (4591)
FB36's recipes (148)

Required Modules

Other Information and Tasks