Fractals/color julia

Visualising complex functions

Color space in HSV-projection

Suppose we want to visualise a complex valued function like

${\begin{aligned}f:\mathbb {C} &\;\to \mathbb {C} \\z&\;\mapsto w=f(z)\\\end{aligned}}$

In order to color $w$ we decompose it into its absolute value $|w|$ and its argument $\arg(w)$.

Then, we assign the color

$\mathrm {HSV} \left({\tfrac {1}{2\pi }}\arg(w),\;1-g_{a_{s},b_{s}}(|w|),\;g_{a_{v},b_{v}}(|w|)\right)$

to the point representing $z$. In this HSV color space, all values are in the range from 0 to  1. The first component (hue) of the HSV color depends only on the argument of $w$ and the second and third component (saturation and value) depend only on the absolute value of $w$. We use a transformation $g_{a,b}$ on $w$ to map it into the interval $[0,1]$. For $g$ see section helper functions.

Examples

The values indexed s (saturation) control the transition saturated colors→white (resp. gray scale), i.e. intermediate values → infinity. The values indexed h (value/valenz) control the transition black→bright, i.e. zero→nonzero. Parameter a controls where the transition takes place: a is just the radius of the circle dividing the two regions (dark/bright, saturated/gray, etc.) Parameter b controls how sharp the transition is: b small = soft, b large = sharp.

The following images all show the range [-10,10]×[-10,10]$i$ and use $a_{s}=5$ (radius of rainbow) and $a_{v}=1$ (radius of black disc).

In the image with swapped meanings of S and V zero is printed in white and infinity in black.

Gallery

• Some examples of colorings
• 
  $b_{s}=b_{v}=10$
• 
  $b_{s}=b_{v}=1$
• 
  $b_{s}=0.05,b_{v}=1/3$
• 
  $b_{s}=0.05,b_{v}=1/3$
  Roles of S and V interchanged

• Analysis of critical point of z→z²+c
• 
  n=1: starting values (color map)
• 
  n=2
• 
  n=3
• 
  n=4
• 
  n=17
• 
  n=18
• 
  n=19
• 
  n=20

Visualising Julia sets

Julia sets and Fatou sets can be nice. Their computergraphical generation can be hard.

In the remainder of this page, we work out a method which can be used to visualise the Julia set of a complex function ƒ. The advantage is that you don't need to know attractors of the iteration

$z\mapsto f(z)$

The generated images will be smooth in the Fatou part.

Overview: Imaging methods

Note that there exist other approaches to color complex dynamical systems like

Inverse Iteration Method (IIM)

Compute the preimages of ƒ, i.e. compute the reverse orbit. Because the stability of the fixed points turns from attractive to repelling and vice versa, one just choses a complex number and looks where it goes under the invers iteration. The trouble is that the results will not uniformly distributet in $J$ and that you have to compute the inverse of ƒ.

Coloring by speed of attraction (CSA)

Taking a value from the lattice of points to color, perform iterations until the iterative is close to an attractor. Color the point according to the number of iterations needed to bring it close enough to the attractor.

This method is commonly used to visualize Julia sets of polynomials and Julia sets that are attached to Newton's famous method for finding the zeros of a function. Polynomials or degree > 1 always have infinity as a super-attractive fixed point. The rational function that occurs in Newton's method has always the function's roots as attractive fixed points. However, in both cases there may be other attractors, which – moreover – need not to consist of just one point.

Escape Time Algorithm (ETA)

If ∞ is an attractor, i.e. a fixed point of the process, then color the point according to the number of iterations – the time – it takes until one sees the point escapes towards ∞. If the point does not escape during the maximum number of iterations, the point is colored as belonging to the Julia set or to the basin of some other attractor. This method works for polynomials. The most prominent Julia sets are the ones for z→z²+c where c is an element of the Mandelbrot set or not far away from the mandelbrot set. If you see a picture of such a Julia set, it is likely that ETA had been used to get the picture.

Cauchy Convergence Algorithm (CCA)

In the remainder of this page, I will present a different approach whose idea as basically the same es that of Escape Time Algorithm. However, no basin of attraction must be known in advance and the different basins of attraction can be separated and be colored differently. The approch uses the notion of Cauchys convergence. Instead of observing the orbit of a point, this method observes how the distance of two nearby points z and z+ε evolves as these two values are iterated. If the difference tends to 0, then the point heads for an attractor. If the difference does not approach 0, then the point is close to (or a part of) the Julia set.

The Metric

Let $\varrho$ be a canonical projection of the compactified complex plane onto the Riemann sphere S₂:

$\varrho :{\overline {\mathbb {C} }}=\mathbb {C} \cup \{\infty \}\rightarrow S_{2}$

This gives us a metric d: As distance between two points in the complex plane we take their distance on the sphere, i.e. the length of the orthodrome. This means that the metric is bounded by π and even the distance to ∞ (which is now the north pole) is finite.

In order to compute the distance between two points z and w we rotate the sphere S₂ in such a way that

1. w maps to 0
2. z maps to the positive real axis

After these transformations the distance can be computed quite easily. The rotation can be accomplished by a squence of isometric Möbius transformations. All an all, we get

$t_{w}:z\mapsto {\begin{cases}|z|\;,&{\text{ if }}w=0\\{\textstyle {\frac {1}{|z|}}}\;,&{\text{ if }}w=\infty \\{\textstyle {\frac {1}{|w|}}}\;,&{\text{ if }}z=\infty \\&\\{\textstyle {\frac {1}{|w|}}}\cdot \left|{\frac {z{\overline {w}}-|w|^{2}}{z{\overline {w}}+1}}\right|\;,&{\text{ else}}\end{cases}}$

which is left as an exercise to the reader. The bar denotes complex conjugation. The metric is then

$d(z,w)=\,\!2\arctan(t_{w}(z))$

The nice feature that is introduced by d is that sequences that formerly diverged against infinity now converge towards infinity, i.e. towards the north pole of S₂.

Stability of Orbits

Recall the definition of the Julia set for a contraction mapping ƒ. The definition implies some facts on the stability of the iteration

$z\mapsto f(z)\mapsto f^{2}(z)\mapsto f^{3}(z)\mapsto \cdots$

where ƒⁿ denotes the n-th iterative of ƒ:

$f^{\,n}=f\circ f\circ \cdots \circ f$

The set

$\{f^{\,n}(z)\}_{n\in \mathbb {N} _{0}}$

is called Orbit of z (under ƒ).

The orbits of two points z and w behave similar − in some sense − if z and w lie close together and belong to the Fatou set F_ƒ which is the complement of the Julia set J_ƒ of ƒ. If z is an element of J_ƒ then z and w will behave quite different, even if w itself is an element of J_ƒ.

To get a notion of the stability of an orbit, we set

$\Sigma _{n}(z)=\Sigma _{n}(z,\varepsilon )=\sum _{k\leqslant n}d\,(f^{k}(z),f^{k}(z+\varepsilon ))$

for a small ε and with the metric d from above. This means that we take two points which are close together, and then we summarize their distances as ƒ makes them jump around on the Riemann sphere. Note that for any fixed ε the sum can diverge for large n even if z is a Fatou point.

However, we can use Σ_n to measure how close a point is to J_ƒ: the larger the sum is, the more instable is the iteration and the closer is the point to the Julia set of ƒ.

To destinguish points of (or close to) the Julia set from points in the Fatou set, we need a creterion. To get it, we compute all the Σ-values for the points that we want to color, i.e. for all points in a lattice Γ. After computing these values, we do a little bit of statistics to get the expected value E and the standard deviation σ for the set of Σ-values: Let

${\begin{aligned}E_{n}&=E_{n}(\Gamma ,\varepsilon )=E\left\{\log(\delta +\Sigma _{n}(z,\varepsilon ))\right\}_{\!z\,\in \,\Gamma }\\\sigma _{n}&={\sqrt {D_{n}}}\\\end{aligned}}$

Remind that

${\begin{aligned}EX&={\textstyle {\frac {1}{|X|}}}\sum \limits _{x\in X}x\\D\!X&=EX^{2}-(EX)^{2}\\\end{aligned}}$

It turns out that the values Σ are widely spread over several scales. Therefore, we do not use Σ directly. Instead, we use the logarithm of Σ. The value δ is just a small constant to avoid the logarithm's input to be zero.

Coloring

Now, we have all we need to color a point:

1. choose
   1. a lattice $\Gamma$ of points $z$ to color
   2. the number $n$ of iterations to perform
   3. a small $\varepsilon$
2. for all points, compute $\Sigma _{n}(z,\varepsilon )$
3. compute $E_{n}$ and $\sigma _{n}$
4. compute

$J(z)=\,\ln(\delta +\Sigma _{n}(z))-K$

for some constant $K$. Because $K$ will be used to separate points that belong to the Julia set ($J(z)>0$) from points that to not ($J<0$), reasonable values for $K$ are greater than $E_{n}$. Try settings like $K=E_{n}+2\sigma _{n}$ or $E_{n}+\sigma _{n}$ or the like. If $J(z)>0$ then we color $z$ as belonging to the Julia set. If $J(z)<0$ we can use that value to shade the Fatou set. If we know some attractor, we can check if ƒ^$n$$(z)$ is close to it and use that information, too.

To map values to the valid ranges for saturation and brightness we use the function $h$ from section helper function h.

Modifications

The computation of $E_{n}$ takes a lot of time. The visualisation process needs two passes:

1. compute $E_{n}$ from all $\Sigma _{n}(z)$
2. color the points using $E_{n}$ and recomputed $\Sigma _{n}(z)$

Alternatively

1. compute and store all $\Sigma _{n}(z)$
2. compute $E_{n}$ and color the points

An other approach looks like that:

Find the smallest $n$ so that $d_{k}(z,z+\varepsilon )$ is below a given bound for all $k\geq n$. If no such $n$ can be found, then assume $z$ to be an element of the Julia set. Otherwise, use $n$ to color the Fatou set.

This algorithm is a variant of the escape time algorithm (ETA). Note that in ETA the point does not really escape (at least if we are on the sphere), it just converges to ∞. This approach is similar. However, we don't need to know an attractive fixed point of ƒ.

Up to now, I didn't try the modified version. One disadvantage may be that the Fatou set will no more appear smooth colored. Then I am not sure if this modification is really an advantage, because the iteration must be done until a given maximum number of iterations is reached. Note that even if $d_{n}$ is under the bound for some $n$ the distance $d_{k}$ can rise again. I cannot say if this effect is crucial or can be neglected...

Gallery

• 
  $z\mapsto z^{2}$
• 
  $z\mapsto z^{2}+i$
• 
  $z\mapsto z^{2}-1$
• 
  ${}_{z\mapsto z^{2}-0.742+i/10}$
• 
  $z\mapsto N_{z^{3}-1}$
• 
  $z\mapsto N_{z^{3}-2z+2}$

$N$ denotes the Newton operator

$N:f\mapsto {\mathit {id}}-{\frac {f}{f'}}$

Using Critical Points Absorption (CPA)

The previous method yeilds neat, smooth colorings and requires least knowledge about the dynamics of the process. However, it is quite time consuming. Teh following approach is an extension of escape time algorithm (ETA) for polynomials.

Let ƒ be a polynomial of degree d > 1 over C. Such a polynomial has at most d critical points: infinity and the at most d–1 zeroes of ƒ′. It is well known that each attractor of the process z→ƒ(z) absorbs at least one critical point. Suppose z_K is a critical value, then ƒⁿ(z_K) comes arbitrarirly close to one of the attractive cycles of ƒ.

A process for a quadratic polynomial ƒ(z) = z² + c is the simplest case: The critical values are 0 and  ∞. As ∞ absorbs itself, only 0 is left, and we have the following cases. M denotes the Mandelbrot set.

$c\notin M$

Easy case. 0 is absorbed by ∞, and J(ƒ) consists just of Cantor dust. Escape time to ∞ can be used to color all points.

$c\in M\setminus \partial M$

If c is an element of the interior of M, then 0 will be absorbed by a (super) attractive cycle. Compute

$w=w(c,n)=f^{n}(0)$

for a sufficiently large n. As w (basically) only depends on c, it can be precomputed before starting the visualization of J(ƒ). Notice that w is the element of a cycle that might have more than one element, i.e. w is only unique modulo that cycle.

Coloring for a point z in C is then as easy as ETA:

• If z is absorbed by ∞, use escape time to color z.
• If ƒⁿ(z) comes close to w, i.e if |ƒⁿ(z) – w| < ε for the first time, use that n to color z.
• If ƒⁿ(z) neither comes close to ∞ nor comes close to w, then color z as part of J(ƒ). Only few z will fall into this category.

Helper functions

Helper function t

some $t\,$

is a smooth, monotone sigmoidal transition from −1 to 1 that satisfies $t'(0)=1$ and $t(x)=-t(-x)$. There are many choices for it. Some of them are

${\begin{aligned}t:\mathbb {R} &\to (-1,1)\\x&\mapsto \tanh \,x\\&\mapsto {\tfrac {2}{\pi }}\arctan \left({\tfrac {\pi }{2}}x\right)\\&\mapsto {\tfrac {2}{\pi }}\operatorname {gd} \left({\tfrac {\pi }{2}}x\right)\\&\mapsto {\frac {x}{\sqrt {1+x^{2}}}}\\&\mapsto {\frac {x}{1+|x|}}\\\end{aligned}}$

with gd denoting the gudermannian function.

• All functions in a Desmos plot.

Helper function h

$h_{a}\,$

Helper function $h_{a}$ is almost the same like $t$ but it maps to $(0,1)$ and we can adjust the speed of transition by parameter $a$.

${\begin{aligned}h_{a}:\mathbb {R} &\to (0,1)\\x\;&\mapsto {\tfrac {1}{2}}+{\tfrac {1}{2}}t(2ax)\end{aligned}}$

with $a>0$. Then $h_{a}$ is symmetric to the point (0, 1/2) and

${\begin{alignedat}{2}h_{a}&(0)&&={\tfrac {1}{2}}\\h'_{a}&(0)&&=a\\h_{a}&(-\infty )&&=0^{+}\\h_{a}&(\infty )&&=1^{-}\\\end{alignedat}}$

Negative values are mapped to values between 0 and 1/2. Positive values are mapped to values between 1/2 and 1. The parameter $a$ controls how fast the transition will be.

If we want a falling function, we can use the symmetry

$h_{a}(x)\,=1-h_{a}(-x)=1-h_{-a}(x)$

i.e. we negate $x$ or $a$.

Helper function g

$g_{a,b}\,$

This function maps the positive numbers to the interval $[0,1)$.

${\begin{aligned}g_{a,b}:\mathbb {R} ^{\geqslant \,0}&\to (0,1)\\x\;&\mapsto {\tfrac {1}{2}}+{\tfrac {1}{2}}\,t(2b\cdot a\cdot u(x/a))\end{aligned}}$

for some function $u$ that is defined below. If $u$ is appropriately chosen then for $g$ the following holds

${\begin{array}{lll}g_{a,b}(0)&=&0\\g_{a,b}(\infty )&=&1\\g_{a,b}(a)&=&1/2\\g'_{a,b}(a)&=&b\\g'_{a,b}(x)&>&0{\text{ for }}x>0\end{array}}$

This means that $g_{a,b}$

• grows continuously from 0 to 1 as $x$ grows from $0$ to $\infty$
• we can control where $g$ crosses the middle between 0 and 1 by specifying parameter $a$
• we can control how fast $g$ passes the point $(a,1/2)$ by specifying parameter $b$

We are left with determining the finction $u$ on $\mathbb {R} ^{>0}$ with

$u\,$

$u$ must satisty

${\begin{array}{lll}u(0^{+})&=&-\infty \\u(\infty )&=&\infty \\u(1)&=&0\\u'(1)&=&1\end{array}}$

For $u$ we set

$u(x)={\tfrac {1}{2}}\left(x-{\tfrac {1}{x}}\right)$

Helper function w

$w_{a}\,$

This function maps the positive numbers to the interval $[0,1)$.

${\begin{aligned}w_{a}:\mathbb {R} ^{\geqslant \,0}&\to [0,1)\\x\;&\mapsto t(ax)\end{aligned}}$

By $a$ we can control its slope in the origin:

$w'(0)=\,a$

• Interactive Desmos plot of $w$

Circular arc through (0,0) and (1,1)

A circular arc through (0,0) and (1,1) that has given slope of $a$ in (0,0) and $-1/a$ in (1,1):

$\operatorname {arc} (x)={\begin{cases}x&;{\text{ if }}a=1\\y_{0}+\operatorname {sign} (\ln a){\sqrt {r^{2}-(x-x_{0})^{2}}}&;{\text{ if }}a\neq 1\\\end{cases}}$

where

${\begin{aligned}x_{0}&={\frac {a}{a-1}}\\y_{0}&=1-x_{0}\\r&={\sqrt {x_{0}^{2}+y_{0}^{2}}}\\\end{aligned}}$

The center of the circle is incident on the line $y=1-x$.

• Desmos plot

See also

• coloring the parameter plane and the Mandelbrot set
• general coloring algorithms: ( feature of the pixel -> color ): $T:scalarvalue\to colorvalue$