Skip to content

Instantly share code, notes, and snippets.

@maltefl
Forked from smukkejohan/chaossearh.py
Created February 21, 2012 18:55
Show Gist options
  • Save maltefl/1878126 to your computer and use it in GitHub Desktop.
Save maltefl/1878126 to your computer and use it in GitHub Desktop.
Python port of Paul Bourke's random strange attractor generator using lyapunov exponents
# Python port of Paul Bourke's http://local.wasp.uwa.edu.au/~pbourke/fractals/lyapunov/gen.c
# By Johan Bichel Lindegaard - http://johan.cc
import math
import random
from PIL import Image, ImageDraw
import argparse
import os
parser = argparse.ArgumentParser(description='Search for chaos.')
#parser.add_argument('-i', dest='maxiterations' metavar='N', type=int,
# help='Maximum iterations.')
args = parser.parse_args()
MAXITERATIONS = 100000
NEXAMPLES = 1000
def createAttractor():
for n in range(NEXAMPLES):
lyapunov = 0
xmin= 1e32
xmax=-1e32
ymin= 1e32
ymax=-1e32
ax, ay, x, y = [], [], [], []
# Initialize coefficients for this attractor
for i in range(6):
ax.append(random.uniform(-2, 2))
ay.append(random.uniform(-2, 2))
# Calculate the attractor
drawit = True;
x.append(random.uniform(-0.5, 0.5))
y.append(random.uniform(-0.5, 0.5))
d0 = -1
while d0 <= 0:
xe = x[0] + random.uniform(-0.5, 0.5) / 1000.0
ye = y[0] + random.uniform(-0.5, 0.5) / 1000.0
dx = x[0] - xe
dy = y[0] - ye
d0 = math.sqrt(dx * dx + dy * dy)
for i in range(MAXITERATIONS):
# Calculate next term
x.append(ax[0] + ax[1]*x[i-1] + ax[2]*x[i-1]*x[i-1] + ax[3]*x[i-1]*y[i-1] + ax[4]*y[i-1] + ax[5]*y[i-1]*y[i-1])
y.append(ay[0] + ay[1]*x[i-1] + ay[2]*x[i-1]*x[i-1] + ay[3]*x[i-1]*y[i-1] + ay[4]*y[i-1] + ay[5]*y[i-1]*y[i-1])
xenew = ax[0] + ax[1]*xe + ax[2]*xe*xe + ax[3]*xe*ye + ax[4]*ye + ax[5]*ye*ye
yenew = ay[0] + ay[1]*xe + ay[2]*xe*xe + ay[3]*xe*ye + ay[4]*ye + ay[5]*ye*ye
# Update the bounds
xmin = min(xmin,x[i])
ymin = min(ymin,y[i])
xmax = max(xmax,x[i])
ymax = max(ymax,y[i])
# Does the series tend to infinity
if xmin < -1e10 or ymin < -1e10 or xmax > 1e10 or ymax > 1e10:
drawit = False
print "infinite attractor"
break
# Does the series tend to a point
dx = x[i] - x[i-1]
dy = y[i] - y[i-1]
if abs(dx) < 1e-10 and abs(dy) < 1e-10:
drawit = False
print "point attractor"
break
# Calculate the lyapunov exponents
if i > 1000:
dx = x[i] - xenew
dy = y[i] - yenew
dd = math.sqrt(dx * dx + dy * dy)
lyapunov += math.log(math.fabs(dd / d0))
xe = x[i] + d0 * dx / dd
ye = y[i] + d0 * dy / dd
# Classify the series according to lyapunov
if drawit:
if abs(lyapunov) < 10:
print "neutrally stable"
drawit = False
elif lyapunov < 0:
print "periodic {} ".format(lyapunov)
drawit = False
else:
print "chaotic {} ".format(lyapunov)
# Save the image
if drawit:
saveAttractor(n,ax,ay,xmin,xmax,ymin,ymax,x,y)
def saveAttractor(n,a,b,xmin,xmax,ymin,ymax,x,y):
width, height = 500, 500
if not os.path.exists("output"):
os.makedirs("output")
# Save the parameters
with open("output/{}.txt".format(n), "w") as f:
f.write("{} {} {} {}\n".format(xmin, ymin, xmax, ymax))
for i in range(6):
f.write("{} {}\n".format(a[i], b[i]))
f.close()
# Save the image
image = Image.new("RGBA", (width, height))
draw = ImageDraw.Draw(image)
for i in range(MAXITERATIONS):
ix = width * (x[i] - xmin) / (xmax - xmin)
iy = height * (y[i] - ymin) / (ymax - ymin)
if i > 100:
draw.point([ix, iy], fill="black")
image.save("output/{}.png".format(n), "PNG")
print "saved attractor to ./output/{}.png".format(n)
createAttractor()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment