Lyapunov Fractal

Lyapunov fractals are mathematical objects living in the plane, graphing regions of stability or chaos regarding a certain logistical map following an ever-changing population growth, alternating between two values.
I implemented a Lyapunov fractal viewer in Java 1.8 which lets you freely move around the plane, set sequences and iteration depths to explore the fractal. Source code can be seen below or downloaded, as can the Java executable (Lyapunov.java, Lyapunov.jar).

Lyapunov[3840x2160](N=100,Zre=3.0,Zim=3.7,Zom=0.3,Seq=AAAAAABBBBBB) — Zircon Zity

Articles on the topic of Lyapunov fractals are a post by Earl Glynn from the year 2000 in which they talk about their Lyapunov fractal generator written in Pascal — a rewrite of a software written in 1992. Also worth reading is A. Dewdney’s article about Mario Markus’ work (Scientific American, September 1991).

My first encounter with this fractal was whilst browsing Wikipedia and stumbling across its Wikipedia entry. Since it was a fractal and looked fairly intricate and thus explorable, I set out to write my own renderer — I chose to implement it in Java, as it is a compiled language with a nice set of GUI libraries. Since Wikipedia listed an algorithm for generating Lyapunov fractal images, I thought an implementation would be fairly straight-forward. However, als always, the devil is in the detail and I quickly noticed that while the short Wikipedia entry looks convincing and well-written at first glance, lacks to answer deeper questions regarding the very topic it is about.

Lyapunov[3840x2160](N=100,Zre=2.294143894223694,Zim=3.4124391882630345,Zom=0.45724737082761846,Seq=ABBABBB) — Eerie Tendrils

The following is a description of the fractal generation algorithm. It is a modified version of the algorithm detailed by Wikipedia which addresses certain ambiguities the original had (more on those later).

Take as input a complex point $(a,b) \in \mathbb{R}^2$ and a sequence $S^*$ consisting of the letters $\text{A}$ and $\text{B}$ ; for example $S^* = \text{AAAAAABBBBBB}$ .
Construct a function $S \colon \mathbb{N}_0 \to \{\text{A},\text{B}\}, n \mapsto S^*_{(n \mod |S^*|)}$ which returns the zero indexed $n$ -th sequence entry and wraps around to the sequence’s start once the sequence’s end is hit.
Construct a function $r \colon \mathbb{N}_0 \to \{a,b\}, n \mapsto \begin{cases}a&\text{if } S(n) = \text{A}\\b&\text{if } S(n) = \text{B}\end{cases}$ which selects either $a$ or $b$ .
Let $x_0=0.5$ and define $x_n=r(n-1)\cdot x_{n-1}\cdot (1-x_{n-1})$ .
Calculate the Lyapunov exponent $\lambda$ as follows.
$\lambda = \lim\limits_{N \to \infty} \frac{1}{N} \cdot \sum\limits_{n=1}^{N} \log_2{|r(n)\cdot(1-2\cdot x_n)|}$
Since calculating the limit as $N$ goes to infinity is rather difficult in practice, one can simply use a sufficiently large $N$ and pretend to have reached infinity.
Output a color value according to $\lambda$ . I chose green for $\lambda < 0$ and blue for $\lambda \geq 0$ , using arbitrarily defined value ranges to map $\lambda$ to an actual color.

Lyapunov[3840x2160](N=1000,Zre=3.8362800000821116,Zim=3.841203134481926,Zom=0.026588814358957543,Seq=AAAAAABBBBBB) — A Spec of Fractal

The following is a Java snippet which implements the algorithm described above.

// choose a or b according to Seq and n
static double r(String Seq, int n, double a, double b) {
	if (Seq.charAt(n%Seq.length()) == 'A') return a; else return b;
}
 
// calculate a pixel color (0x00rrggbb) for given parameters
static int LyapunovPixel(double a, double b, String Seq, int N) {
	// array of all x_n; x_0, the starting value, initialize all x_n values
	double[] X = new double[N+1]; X[0] = .5;
	for (int n = 1; n <= N; n++) X[n] = r(Seq,n-1,a,b)*X[n-1]*(1-X[n-1]);
	
	// calculate the Lyapunov exponent (to a certain precision dictated by N)
	double lmb = 0;
	for (int n = 1; n <= N; n++)
		lmb += Math.log(Math.abs(r(Seq,n,a,b)*(1-2*X[n])))/Math.log(2);
	lmb /= N;

	// infinity was hit, use a black pixel
	if (Double.isInfinite(lmb)) return 0x000000;

	// color pixel according to Lyapunov exponent
	double MIN = -1, MAX = 2;
	if (lmb < 0) return ((int)(lmb/MIN*255))<<8;
	else         return ((int)(lmb/MAX*255))<<0;
}

Lyapunov[3840x2160](N=100,Zre=2.2222222222222223,Zim=3.0,Zom=1.0,Seq=AB) — Lyapunov Spike

Coming back to Wikipedia’s algorithm, there were a few things I found irritating at best when attempting my implementation and thus addressed in the algorithm description seen above. A closer look at potentially misleading or ambiguos statements follows, together with my reasoning to resolve them.

It is not clear whether the sequence is zero or one indexed, though it has to be zero indexed as $x_0=0.5, x_{n+1}=r_n\cdot\dots$ evaluates to $x_1=r_0\cdot\dots$ , implying the definition of $r_0$ and thus $S_0$ .

It is not clear which logarithm is meant; the natural logarithm, a generic $\log_b$ or the dyadic logarithm — the latter is actually used. To find the actual logarithm base, I dug deeper and used Wikipedia’s external links to find a post by Earl Glynn, answering this question.

It is not clear what is meant by the second half of the sentence beneath the Lyapunov exponent. It reads “… dropping the first summand as $r_0(1-2x_0) = r_n\cdot 0 = 0$ for $x_0 = 0.5$ .” As the sum starts with $n=1$ and one would not dare to calculate $\log_2{|r(0)\cdot(1-2\cdot x_0)|} = \log_2{0}$ for $x_0 = 0.5$ , this sentence’s existence bewilders me.

It is not clear how exactly the colors are derived, only that they have something to do with the Lyapunov exponent. I simply chose arbitrary values that looked convincing for mapping $\lambda$ to a pixel color.

Lyapunov[3840x2160](N=100,Zre=3.015244151652456,Zim=3.692602005918367,Zom=0.05559060566555525,Seq=AAAAAABBBBBB) — Dark Swirl

One of the most frustrating bugs I came across was an unexplainable axis flip. My code generated the fractal just fine except for the fact that every image was flipped along the diagonal crossing the origin with a $45^\circ$ angle to the horizontal. It was as though the coordinates $(a,b)$ were swapped somewhere in my code and I simply could not figure out where.

Finally, after hours of looking at the same code over and over again, a closer look at Earl Glynn’s post brought an end to my misery. Just below the three welcoming fractal renderings, a screenshot of their software is shown — complete with a blue and red line indicating the coordinate system’s orientation. $a$ and $b$ are — contrary to all coordinate systems involving parameters named after the first two letters in the latin alphabet — indeed flipped. Wikipedia’s images must have simply ran with this decision, without noting it.
Because of this flip, when one wants to render images specified in the reversed sequence format, they simply have to swap all letters (for example $\text{BBBABA}$ becomes $\text{AAABAB}$ ).

As Glynn themself says, “I would have reversed the ‘a’ and ‘b’ labels to be consistent with normal ‘x’ and ‘y’ axes conventions, but conform here with the same convention as used by Markus.” (Referring to Mario Markus, co-author of Lyapunov Exponents of the Logistic Map with Periodic Forcing.)

Lyapunov[3840x2160](N=10,Zre=-1.595906489339523,Zim=-1.1504753527142757,Zom=0.5233476330273611,Seq=ABBABBB) — Slurping Cell

After having eliminated all vagueness regarding the fractal generation, writing the outer Java viewer hull was a fairly easy task. As a template I took my Mandelbrot Set III viewer (the reason why most variable names reference complex numbers) complete with multithreading, allowing a pixely fractal exploration until one stops their exploration, letting the thread catch up and display a higher resolution image. The same is done for rendering 4K images and saving them to files in the background — 4K renderings are saved as .png files and named according to their parameters. A list of program controls follows.

Mouse dragging lets one pan the complex plane.
Mouse scrolling lets one zoom the complex plane.
Pressing N evokes a dialogue box where one can specify an iteration depth $N$ .
Pressing S evokes a dialogue box where one can enter a sequence $S^*$ .
Pressing R resets all fractal parameters to the default.
Pressing P initiates a 4K fractal rendering. The 4K fractal rendering thread can be killed by exiting the program prematurely, thus losing (!) the rendering.

// Java 1.8 code; written by Jonathan Frech, 2017

// === About ===
// This Java program renders the Lyapunov fractal.
// Blog post: https://jonathanfrech.wordpress.com/2017/12/30/lyapunov-fractal/

// === Edit History ===
// Fractal viewer core: 22nd, 23rd of July, 23rd, 24th of December 2017
// Lyapunov fractal: 27th, 29th, 30th of December 2017

// === Compilation ===
// $ javac Lyapunov.java; jar -cvfe Lyapunov.jar Lyapunov *.class

// import
import javax.swing.*;
import java.awt.*;
import java.awt.event.*;
import java.awt.image.BufferedImage;
import java.awt.image.DataBufferInt;
import java.awt.Color;
import java.lang.Math.*;
import java.io.File;
import java.io.IOException;
import javax.imageio.ImageIO;
import java.util.ArrayList;
import java.util.Arrays;

class FractalRenderer implements Runnable {
	private MainPanel parent;
	private Thread thread;

	FractalRenderer(MainPanel p) { parent = p; thread = new Thread(this, "FractalRenderer"); thread.start(); }

	public void run() {
		// loop indefinitely, check if fractal has to be drawn
		while (true) {
			// parent dirty, rendering needs to be done
			if (parent.dirty) {
				// save parameters, render, make sure that parameters have not changed while rendering
				int N = parent.N; double Zre = parent.Zre, Zim = parent.Zim, Zom = parent.Zom; String Seq = parent.Seq;
				LyapunovRenderer.render(parent.WIDTH, parent.HEIGHT, parent.pxl, N, Zre, Zim, Zom, Seq);
				if (N == parent.N && Zre == parent.Zre && Zim == parent.Zim && Zom == parent.Zom && Seq.equals(parent.Seq)) { parent.dirty = false; parent.repaint(); }
			}

			// did not need to render, sleep
			else { try { Thread.sleep(100); } catch (InterruptedException e) {} }
		}
	}
}

// store fractal parameters
class FractalParameters {
	int N; double Zre, Zim, Zom; String Seq;
	FractalParameters(int N, double Zre, double Zim, double Zom, String Seq) { this.N = N; this.Zre = Zre; this.Zim = Zim; this.Zom = Zom; this.Seq = Seq; }
}

class FractalRendererQueue implements Runnable {
	private Thread thread;
	private final int WIDTH, HEIGHT;
	private BufferedImage img;
	private int[] pxl;
	private ArrayList<FractalParameters> queue = new ArrayList<>();

	FractalRendererQueue(int W, int H) {
		WIDTH = W; HEIGHT = H;
		img = new BufferedImage(WIDTH, HEIGHT, BufferedImage.TYPE_INT_RGB);
		pxl = ((DataBufferInt) img.getRaster().getDataBuffer()).getData();
		thread = new Thread(this, "FractalRendererQueue");
		thread.start();
	}

	// add another item to the queue
	void add(FractalParameters fp) { queue.add(fp); }

	public void run() {
		// loop indefinitely, check if fractal has to be drawn
		while (true) {
			// item in queue, render
			if (queue.size() > 0) {
				FractalParameters fp = queue.remove(0);
				LyapunovRenderer.render(WIDTH, HEIGHT, pxl, fp);
				try {
					String fn = System.getProperty("user.dir") + "/Lyapunov["+WIDTH+"x"+HEIGHT+"](N="+fp.N+",Zre="+fp.Zre+",Zim="+fp.Zim+",Zom="+fp.Zom+",Seq="+fp.Seq+").png";
					ImageIO.write(img, "png", new File(fn));
				} catch (Exception e) { }
			}

			// queue empty, sleep
			else { try { Thread.sleep(100); } catch (InterruptedException e) {} }
		}
	}
}

class LyapunovRenderer {
	// r_n
	static double r(String Seq, int n, double a, double b) {
		if (Seq.charAt(n%Seq.length()) == 'A') return a;
		else return b;
	}

	// render a Lyapunov fractal on given pixel array
	static void render(int W, int H, int[] pxl, int N, double Zre, double Zim, double Zom, String Seq) {
		// array of all x_n; x_0, the starting value
		double[] X = new double[N+1]; X[0] = .5;

		// complex height and width, ln(2) constant
		double h = Zom*2, w = h*W/H;
		double ln2 = Math.log(2);

		// used for coloring
		double MIN = -1, MAX = 2;

		// loop through pixels
		for (int y = 0; y < H; y++) { for (int x = 0; x < W; x++) {
			// complex coordinates
			double a = Zre-w/2.+w*x/W, b = Zim+h/2.-h*y/H;

			// only focus on a certain subset of the complex plane (disabled)
			//if (a < 0 || a > 4 || b < 0 || b > 4) { pxl[x+y*W] = 0; continue; }

			// initialize all x_n values
			for (int n = 1; n <= N; n++) X[n] = r(Seq, n-1, a, b)*X[n-1]*(1-X[n-1]);

			// calculate the Lyapunov exponent (to a certain precision dictated by N)
			double lmb = 0;
			for (int n = 1; n <= N; n++) lmb += Math.log(Math.abs(r(Seq, n, a, b)*(1-2*X[n])))/ln2;
			lmb /= N;

			// infinity was hit, use a black pixel
			if (Double.isInfinite(lmb)) { pxl[x+y*W] = 0; continue; }

			// color pixel according to Lyapunov exponent and defined MIN, MAX values (colors hard coded)
			if (lmb < 0) pxl[x+y*W] = ((int)(lmb/MIN*255))<<8;
			else         pxl[x+y*W] = ((int)(lmb/MAX*255))<<0;
		}}
	}

	// render using a FractalParameters class to hold fractal parameters
	static void render(int W, int H, int[] pxl, FractalParameters fp) { render(W, H, pxl, fp.N, fp.Zre, fp.Zim, fp.Zom, fp.Seq); }
}

class MainPanel extends JPanel {
	// image dimensions
	final int WIDTH_SMALL = 64, HEIGHT_SMALL = 48;
	final int WIDTH = 640, HEIGHT = 480;
	final int WIDTH_LARGE = 3840, HEIGHT_LARGE = 2160;

	// buffered images and their pixel arrays
	BufferedImage img_small, img;
	int[] pxl_small, pxl;

	// fractal parameters
	boolean dirty = false;
	double Zom, Zre, Zim;
	String Seq; int N;

	// queue thread
	private FractalRendererQueue Renderer;

	// mouse paning
	private boolean paning = false;
	private int paningx, paningy;
	private double paningre, paningim;
	private final double panspeed = 50;

	// zoom change, maximal zoom
	private final double Z = .9;
	private final double ZomMax = 10;

	// constructor
	MainPanel() {
		// initialize images
		img_small = new BufferedImage(WIDTH_SMALL, HEIGHT_SMALL, BufferedImage.TYPE_INT_RGB);
		img = new BufferedImage(WIDTH, HEIGHT, BufferedImage.TYPE_INT_RGB);
		pxl = ((DataBufferInt) img.getRaster().getDataBuffer()).getData();
		pxl_small = ((DataBufferInt) img_small.getRaster().getDataBuffer()).getData();

		// start threads
		new FractalRenderer(this);
		Renderer = new FractalRendererQueue(WIDTH_LARGE, HEIGHT_LARGE);

		// add listeners, focus, reset and thereby start program
		addlisteners(); setFocusable(true); reset();
	}
	
	// reset parameters
	private void reset() {
		// Zircon Zity
		Seq = "AAAAAABBBBBB"; Zom = .3; Zim = 3.7; N = 100;
		double h = Zom*2, w = h*WIDTH/HEIGHT;
		Zre = 3.4-w/2;
		
		draw();
	}

	// draw preview, tell thread to draw
	private void draw() {
		if (Zom > ZomMax) Zom = ZomMax;
		LyapunovRenderer.render(WIDTH_SMALL, HEIGHT_SMALL, pxl_small, N, Zre, Zim, Zom, Seq);
		dirty = true; repaint();
	}

	// draw component
	@Override
	public void paintComponent(Graphics g) {
		super.paintComponent(g);

		if (dirty) g.drawImage(img_small, 0, 0, WIDTH, HEIGHT, null);
		else g.drawImage(img, 0, 0, WIDTH, HEIGHT, null);

		g.dispose();
	}

	// component size
	@Override
	public Dimension getPreferredSize() { return new Dimension(WIDTH, HEIGHT); }

	// add listeners
	private void addlisteners() {
		addKeyListener(new KeyAdapter() {
			public void keyPressed(KeyEvent e) {
				// complex frame size
				double im = Zom/2, re = im*WIDTH/HEIGHT;

				// set N
				if (e.getKeyCode() == KeyEvent.VK_N) {
					try {
						String I = (String) JOptionPane.showInputDialog(null, "N = ", "Set N (currently "+N+")", JOptionPane.QUESTION_MESSAGE, null, null, Integer.toString(N));
						if (N > 0) { N = Integer.parseInt(I); draw(); }
					}
					catch (NullPointerException ne) { }
					catch (NumberFormatException ne) { }
				}

				// set Seq
				else if (e.getKeyCode() == KeyEvent.VK_S) {
					try {
						String I = (String) JOptionPane.showInputDialog(null, "Seq = ", "Set Seq (currently "+Seq+")", JOptionPane.QUESTION_MESSAGE, null, null, Seq);
						String nSeq = "";
						for (char c : I.toCharArray()) {
							if (c == 'A' || c == 'a') nSeq += "A";
							if (c == 'B' || c == 'b') nSeq += "B";
						}
					  if (nSeq.length() > 0) { Seq = nSeq; draw(); }
					} catch (NullPointerException ne) { }
				}

				// reset
				else if (e.getKeyCode() == KeyEvent.VK_R) reset();

				// save fractal rendering
				else if (e.getKeyCode() == KeyEvent.VK_P) { Renderer.add(new FractalParameters(N, Zre, Zim, Zom, Seq)); }
			}
		});

		addMouseListener(new MouseAdapter() {
			public void mousePressed(MouseEvent e) {
				paning = true;
				paningx = e.getX(); paningy = e.getY();
				paningre = Zre; paningim = Zim;
			}

			public void mouseReleased(MouseEvent e) {
				if (paning) {
					if (dirty) pan(e);
					paning = false;
				}
			}
		});

		addMouseMotionListener(new MouseAdapter() {
			public void mouseDragged(MouseEvent e) { pan(e); }
		});

		addMouseWheelListener(new MouseAdapter() {
			public void mouseWheelMoved(MouseWheelEvent e) {
				// scroll intensity, complex frame size
				int n = e.getWheelRotation();
				double im = Zom*2, re = im*WIDTH/HEIGHT;

				// do not scroll (and thereby render) unnecessarily
				if (n > 0 && Zom == ZomMax) return;

				// mouse cursor's real and imaginary part
				double r = Zre-re/2+re*e.getX()/WIDTH;
				double i = Zim-im/2+im*(HEIGHT-e.getY())/HEIGHT;

				// zoom
				Zom *= Math.pow(Z, -n);
				if (Zom > ZomMax) Zom = ZomMax;
				im = Zom*2; re = im*WIDTH/HEIGHT;

				// shift so that cursor is at its previous relative position
				Zre = re/2-re*e.getX()/WIDTH+r;
				Zim = im/2-im*(HEIGHT-e.getY())/HEIGHT+i;

				// redraw
				draw();
			}
		});
	}

	// shift frame based on mouse dragging
	private void pan(MouseEvent e) {
		// calculate pixel delta, complex delta
		int dx = paningx-e.getX(), dy = paningy-e.getY();
		double im = Zom*2, re = im*WIDTH/HEIGHT;
		double dre = re*dx/WIDTH, dim = im*dy/HEIGHT;

		// add complex delta
		Zre = paningre + dre;
		Zim = paningim - dim;

		// redraw
		draw();
	}
}

class Lyapunov extends JFrame {
	private Lyapunov() {
		setTitle("Lyapunov"); setDefaultCloseOperation(WindowConstants.EXIT_ON_CLOSE);
		setLocationByPlatform(true); setResizable(false);
		add(new MainPanel()); pack(); setVisible(true);
	}

	public static void main(String[] args) { new Lyapunov(); }
}

One thought on “Lyapunov Fractal”

Pingback: Third Anniversary – J-Blog

Like this:

One thought on “Lyapunov Fractal”

Leave a Reply Cancel reply