Variable exponent in Mandelbrot algorithm

Animation of Mandelbrot with variable exponent

In my previous two posts, I generated images of the Mandelbrot Set using the conventional iterative mapping function to generate the sequence of complex values for each point in the area of the complex plane that was shown. I’m struggling to develop a better intuitive sense of where the form and complexity of the Mandelbrot set originates – why on Earth is it that crazy shape? I’m particularly curious about the sensitivity of the process on different aspects of the mapping function used. What conditions are necessary to produce a fractal pattern?

So I thought I’d begin by varying the exponent (its conventional value is 2). The animation above shows how the set varies as the value of the exponent in the mapping function is varied from 1.0 to 3.0 in 100 steps.

One thing that really surprised me about the images produced when the exponent is between 1 and 2 is how “edgy” the patterns are. They are reminiscent of fractured tiling in a way that the classic Mandelbrot pattern is not. I’ll include an additional zoomed in colour image at the bottom that shows what I mean.

This is the C program I used to generate each frame of the animation:

//
// mandelbrot.c - Visualise the Mandelbrot Set
// written by Ted Burke - last updated 16-12-2012
//
// The command line arguments are the real and imaginary parts of
// the centre value, the real axis range and the output filename.
// 
// To compile with MinGW:
//
//		gcc mandelbrot.c -o mandelbrot.exe -lm
//
// To run (this example is centred on point x=-0.7485,y=0.1):
//
//		mandelbrot.exe -0.7485 0.1 0.0001 ms.pgm
//

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>

#define PI 3.14159265

#define W 640
#define H 640

// Pixel array
int p[H][W];

int main(int argc, char *argv[])
{
	int x, y;
	double _Complex cmin, crange, c, z;
	double real_range;
	double limit = 10.0;
	double maxval = 0;
	double power;
	
	// Parse command line arguments
	if (argc != 6)
	{
		printf("Usage: mandelbrot LEFT_REAL TOP_IMAG REAL_RANGE POWER FILENAME\n");
		return 0;
	}
	
	cmin = atof(argv[1]) + atof(argv[2])*_Complex_I;
	real_range = atof(argv[3]);
	crange = real_range + (H*real_range/W)*_Complex_I;
	power = atof(argv[4]);
	
	// Calculate pixel values
	for (y=0;y<H;++y)
	{
		for (x=0;x<W;++x)
		{
			// Calculate c value for this pixel
			c = cmin + (x/(double)W)*creal(crange) + (y/(double)H)*cimag(crange)*_Complex_I;
			
			// Iterate for current c value
			p[y][x] = 255;
			z = 0;
			while (cabs(z) < limit && p[y][x] > 0)
			{
				z = cpow(z, power) + c;
				p[y][x] *= 0.98;
			}
			maxval = (p[y][x] > maxval) ? p[y][x] : maxval;
		}
		
		fprintf(stderr, "\rRow %d calculated", y);
	}
	fprintf(stderr, "\r                                   ");
	
	// Output image to PGM file
	FILE *f = fopen(argv[5], "w");
	int r, g, b;
	fprintf(f, "P5\n# Image\n%d %d\n255\n", W, H);
	for (y=0;y<H;++y)
	{
		for (x=0;x<W;++x)
		{
			fputc(p[y][x], f);
		}
		
		fprintf(stderr, "\rRow %d written", y);
	}
	fprintf(stderr, "\n");
	fclose(f);
	
	return 0;
}

This is the Python script I used to create the gif animation:

#
# makegif.py - Generate an animated Mandelbrot gif
# Written by Ted Burke - last updated 16-12-2012
#

import math
import os

x, y, r = -1.5, -1.5, 3.0
p1, p2 = 1.0, 3.0
N = 100

pscale = pow(p2/p1, 1/float(N-1))
p = p1

for n in range(N):
	p = p1 * math.pow(pscale, n)
	
	print "frame %d/%d: " % (n+1, N), x, y, r, p
	
	command = "./mandelbrot " + \
		str(x) + " " + str(y) + " " + \
		str(r) + " " + str(p) + " " + \
		"frame%03d.pgm" % (n)
	os.system(command)
	
	command = "convert frame%03d.pgm frame%03d.png" % (n,n)
	os.system(command)
		
	command = "rm frame%03d.pgm" % (n)
	os.system(command)

os.system("convert -delay 20 frame* ms.gif")

Here’s a zoomed in colour image that illustrates that edgy characteristic of the patterns produced by exponents between 1 and 2.

Zoomed in Mandelbrot image with exponent set to 1.3

I generated that image using this command:

./mandelbrot_colour 0 0.25 0.1 1.3 fractured_colour.ppm

This was the modified C code I used:

//
// mandelbrot_colour.c - Visualise the Mandelbrot Set
// written by Ted Burke - last updated 16-12-2012
//
// The command line arguments are the real and imaginary parts of
// the centre value, the real axis range and the output filename.
// 
// To compile with MinGW:
//
//		gcc mandelbrot.c -o mandelbrot.exe -lm
//
// To run (this example is centred on point x=-0.7485,y=0.1):
//
//		mandelbrot.exe -0.7485 0.1 0.0001 ms.ppm
//

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>

#define PI 3.14159265

#define W 640
#define H 640

// Pixel array
int p[H][W];

int main(int argc, char *argv[])
{
	int x, y;
	double _Complex cmin, crange, c, z;
	double real_range;
	double limit = 10.0;
	double maxval = 0;
	double power;
	
	// Parse command line arguments
	if (argc != 6)
	{
		printf("Usage: mandelbrot LEFT_REAL TOP_IMAG REAL_RANGE POWER FILENAME\n");
		return 0;
	}
	
	cmin = atof(argv[1]) + atof(argv[2])*_Complex_I;
	real_range = atof(argv[3]);
	crange = real_range + (H*real_range/W)*_Complex_I;
	power = atof(argv[4]);
	
	// Calculate pixel values
	for (y=0;y<H;++y)
	{
		for (x=0;x<W;++x)
		{
			// Calculate c value for this pixel
			c = cmin + (x/(double)W)*creal(crange) + (y/(double)H)*cimag(crange)*_Complex_I;
			
			// Iterate for current c value
			p[y][x] = 255;
			z = 0;
			while (cabs(z) < limit && p[y][x] > 0)
			{
				z = cpow(z, power) + c;
				p[y][x] *= 0.98;
			}
			maxval = (p[y][x] > maxval) ? p[y][x] : maxval;
		}
		
		fprintf(stderr, "\rRow %d calculated", y);
	}
	fprintf(stderr, "\r                                   ");
	
	// Output image to PPM file
	FILE *f = fopen(argv[5], "w");
	int r, g, b;
	fprintf(f, "P6\n# Image\n%d %d\n255\n", W, H);
	for (y=0;y<H;++y)
	{
		for (x=0;x<W;++x)
		{
			b = 255.0 * (0.5 + 0.5*cos(2*PI * p[y][x] / maxval));
			r = 255.0 * (0.5 + 0.5*cos(2*PI * p[y][x] / maxval + 2*PI/3.0));
			g = 255.0 * (0.5 + 0.5*cos(2*PI * p[y][x] / maxval + 2*PI/1.5));
			
			//r = p[y][x]; g = p[y][x]; b = p[y][x];
			
			//fprintf(f, "%03d %03d %03d ", r, g, b);
			fputc(r, f);
			fputc(g, f);
			fputc(b, f);
		}
		
		fprintf(stderr, "\rRow %d written", y);
	}
	fprintf(stderr, "\n");
	fclose(f);
	
	return 0;
}
Advertisements
This entry was posted in Uncategorized and tagged , , , , , , , , . Bookmark the permalink.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s