Gesamtschule AachenBrand, Rombachstrasse 99, 52078 Aachen, Germany.
Mandel: software for real and complex dynamics
Mandel is an interactive program for drawing the Mandelbrot set
and Julia sets, and for illustrating and researching their mathematical
properties. It is available on Linux, Windows, and Mac with
a graphical user interface based on the c++ toolkit
Qt
by Digia. Features include:

Draw the Mandelbrot set and Julia sets in two windows at the same time.
Available methods of drawing include escape time, hyperbolic components or
attracting basins, distance estimate, marking external rays or puzzlepieces.
Display Julia sets on the Riemann
sphere.

Trace equipotential lines and external rays, in the parameter plane
and in the dynamic plane as well.

Discuss the combinatorics of external angles, kneading sequences, and
internal addresses. Illustrate the Spider Algorithm and Thurston Algorithm.

Find centers and Misiurewicz points in the parameter plane, or
periodic and preperiodic points in the dynamic plane.

Illustrate renormalization
by marking small Mandelbrot sets and small Julia sets.
Mark
embedded Julia sets
by a kind of preperiodic renormalization.

Display the asymptotic selfsimilarity at
Misiurewicz points
on multiple scales, and the local similarity
to Julia sets.

Draw the parameter plane and Julia sets for other oneparameter
families, which satisfy critical relations or show a
persistent Siegel disk.
Available are unicritical polynomials, BrannerFagella
polynomials, various families of cubic and quartic polynomials,
quadratic rational mappings,
Newton's method,
transcendental functions, a simulation of quasiconformal surgery, the
Tricorn,
Henon mappings, and the Barnsley–Bousch–Thurston IFS.

Compute and display the core entropy and biaccessibility dimension for dyadic angles.

Save images as PostScript *.eps, save or load a *.png image.

Animated demos
are providing an extensive introduction to basic and
advanced aspects of complex dynamics.

Start Mandel from the GAP package
FR by
Laurent Bartholdi.
E.g., this package can compute a mating from combinatorial data.

Future plans include the construction of Hubbard trees and incorporating real iteration.

The program contains a few languages: English, German, Polish
(by Adam Majewski),
and Portuguese
(by Atractor).
The latest release is Mandel 5.10 of January 12, 2014.
The source code is available under the
GNU General Public Licence.
Compile it yourself on Linux, Windows, or Mac after installing
Qt 45.
For your convenience, a Windows executable is offered in addition.
Download the source code (485 kB)
or the Windows exe (5.5 MB).
There is an executable for Mac as well
(7.6 MB).
Download Juliette 1.0 here, a Java applet
providing external rays, renormalization, asymptotic and local similarity.
Other sites with programs on complex dynamics or fractals:
Arnaud Chéritat,
Curt McMullen,
Mitsuhiro Shishikura,
Xaos,
Fractint.
Just for fun ... a PostScript program for computing M.
You may change parameters by editing the source code.
Download mc.eps (1.4 kb).
The DOSprogram logistic.exe deals with the real iteration of z^{2} + c,
which is conjugate to the standard logistic map. It provides features like
drawing qcurves, finding hyperbolic intervals or periodic points, cobweb
graphics. This program may be used and redistributed without restriction.
It comes with no warranty.
Download logistic.exe
2.0 of July 30, 1998. Preliminary manual logistic.ps.
How to draw external rays:
Background:
In the iteration of quadratic polynomials f_{c}(z) = z^{2} + c ,
the parameter c is fixed and the dynamic variable z is evolving when the
mapping is applied again and again. Both c and z are
complex numbers.
The filled Julia set K_{c} contains those values z that are not
going to the
point at ∞
under iteration. The dynamics is chaotic on its boundary ∂K_{c} ,
which is the Julia set.
The mapping f_{c}(z) = z^{2} + c is conjugate to
F(z) = z^{2} in a neighborhood of ∞. The conjugacy is provided
by the Boettcher mapping Φ_{c}(z), which satisfies
Φ_{c}(f_{c}(z)) = F(Φ_{c}(z)). We have
(1) Φ_{c}(z) = lim_{n→∞}
(f_{c}^{n}(z))^{1/2n}
for a suitable choice of the 2^{n}th root. This is made precise by
the following formula, which is valid for large z with the principal value
of the roots:
(2) Φ_{c}(z) = z ∏_{n=1}^{∞}
(1 + c/(f_{c}^{n1}(z))^{2})
^{1/2n}
= z ∏_{n=1}^{∞}
(f_{c}^{n}(z)/(f_{c}^{n}(z)c))
^{1/2n}.
In the dynamics of quadratic polynomials, there is a basic dichotomy:

When the critical orbit stays bounded, then the Julia set is connected, and
Φ_{c}(z) extends to a conformal mapping from its exterior to the
exterior of the unit disk. External rays and equipotential lines in the
dynamic plane are defined as preimages of straight rays and circles,
respectively. By definition, the parameter c belongs to the Mandelbrot set M.

When the critical orbit escapes to ∞, the Julia set is totally
disconnected, and the parameter c is not in M. Now Φ_{c}(z) does
not extend to the critical point z = 0, e.g., but it is welldefined at the
critcal value z = c.
The mapping Φ_{M}(c) is defined by evaluating the Boettcher mapping
at the critcal value z = c, i.e., Φ_{M}(c) = Φ_{c}(c).
Adrien Douady
has shown that it is a conformal mapping from the exterior of the Mandelbrot set
M to the exterior of the unit disk, thus M is connected.
By means of Φ_{M}(c), external rays and equipotential lines are
defined in the parameter plane as well. See the image above. When the argument
θ of an external ray is a rational multiple of 2π, then the ray is
landing at a distinguished boundary point: a rational dynamic ray is landing at
a periodic or preperiodic point in the Julia set. A rational parameter ray is
landing at a Misiurewicz point or a root point in the boundary of the
Mandelbrot set. See the following papers by
Dierk Schleicher
(arXiv:math.DS/9711213) and
John Milnor
(arXiv:math.DS/9905169).
Computation of the external argument:
External rays are characterized by the fact that the argument arg_{c}(z) of
Φ_{c}(z), or the argument arg_{M}(c) of Φ_{M}(c),
equals a given value θ. For large z, formula (2) implies
(3) arg_{c}(z) = arg(z) + ∑_{n=1}^{∞}
1/2^{n} arg(1 + c/(f_{c}^{n1}(z))^{2})
= arg(z) + ∑_{n=1}^{∞}1/2^{n}
arg(f_{c}^{n}(z)/(f_{c}^{n}(z)c)).
Here arg(z) denotes the
argument
of the
complex number
z, i.e., its angle to the positive real axis. This angle is defined only up to
multiples of 2π. E.g., a point z on the negative real axis has the arguments
±π, ±3π, ±5π ... The principal value of the
argument is the unique angle with π<arg(z)≤π. This definition
is used because a function should be singlevalued. However, a discontinuity
is introduced artificially when a point z is crossing the negative real axis,
say going from i through 1 to i: its argument should go from π/2 through
π to 3π/2, but the principal value goes from π/2 to π, jumps to
π, and goes to π/2.
Formula (3) means that the orbit z_{n}=f_{c}^{n}(z) of
z and the arguments arg(z_{n}/(z_{n}c))
are computed. The function arg(z/(zc)) is discontinuous
on the straight line segment from the critical point z = 0 to the critical
value z = c, which is drawn in red in the figure: when z belongs to this line
segment, then z/(zc) will be on the negative real axis. When the Julia set is
connected, the location of the discontinuity can be moved onto a kind of curve
joining 0 and c within the Julia set, such that the argument is continuous in
the whole exterior:

Suppose that z is crossing the red line from right to left, and moving to the
fixed point α_{c} and beyond. The principal value arg(z/(zc))
is jumping from π to π at the red line, goes up to about 0.68π at
α_{c} and stays continuous there.

The modified function arg(z/(zc)) is continuous with value π, as z is
crossing the red line. It goes up to about 1.32π when z is approaching
α_{c} from the right, and jumps to about 0.68π as z is
crossing α_{c} to the left.
When the modified function arg(z/(zc)) is used instead of the principal value
for computing arg_{c}(z) according to (3), this formula will be valid
not only for large z but everywhere outside the Julia set. In practice, the
Julia set may be approximated by two straight line segments from the fixed
point α_{c} to 0 and c, respectively: when z_{n} is in
the triangle between 0, c, and α_{c}, adjust the principal value
of arg(z_{n}/(z_{n}c)) by ±2π. The triangle is given
by three simple
inequalities, cf. mndlbrot::turn(). In the two images below, the
different colors indicate values of arg_{c}(z), and thus, external rays.
The computed value of arg_{c}(z), and thus the color, is jumping on
certain curves joining preimages of 0, where some iterate z_{n} is
crossing the line where arg(z/(zc)) is discontinuous. In the left image, the
principal value of arg(z/(zc)) is used, and in the right image, the argument
is adjusted. The erroneous discontinuities of the computed arg_{c}(z)
are moved closer to the Julia set.
The same modification seems to work for the computation of the argument
arg_{M}(c) in the parameter plane as well, although the Julia set is
disconnected. The left image below shows the discontinuities arising from the
principal value of arg(z/(zc)) in (3), they are appearing on curves joining
centers of different periods. In the right image, the argument is adjusted,
and the discontinuities are moved closer to the Mandelbrot set.
Drawing “all” rays:
To get an overview of all rays, compute the argument at every pixel and map it
to a range of colors. The images above where obtained with Mandel's algorithm
No. 6. The code of mndlbrot::turn() has been
ported to Java by
Evgeny Demidov.
See also the algorithms 7 and 8, where no discontinuity appears, but only the
ends of certain rays are visible.
Adam Majewski
has suggested that boundary scanning (below) can be used to draw several rays simultaneously.
Drawing a single ray:

In the dynamic plane, external rays can be drawn by backwards iteration:
It is most effective for a periodic or preperiodic angle. You must keep track
of points on the finite collection of rays with angles θ, 2θ,
4θ ... Say z_{l,j} corresponds to a radius
R^{1/2l} and the angle 2^{j}θ. Then
f_{c}(z) maps z_{l,j} to z_{(l1),(j+1)}. This point,
which was constructed before, has two preimages under f_{c}(z) . The
one that is closer to z_{(l1),j} is the correct one.
This criterion was proved by
Thierry Bousch.
The ray will look better when you introduce intermediate points. Backwards
iteration is used in Mandel again since version 5.3, both in the plane and
on the sphere.
The current implementation consists of QmnPlane::backRay().
The following methods apply with small modifications to the dynamic plane and
the parameter plane as well:

Choose a sequence of radii tending to 1, e.g.,
R_{n} = R^{1/2n} with R large.
Search corresponding points on the ray, which satisfy
Φ_{c}(z_{n}) = R_{n} e^{iθ} or
Φ_{M}(c_{n}) = R_{n} e^{iθ},
respectively. To this end, the Newton method is employed. Using the
approximation (1), the function N(z) = z  P(z)/P'(z)
is iterated to find a zero of
P(z) = f_{c}^{n}(z)  R e^{i2nθ},
or analogously in the parameter plane. Intermediate points should be used
to make the curve between the points smooth, and to reduce the danger of
the Newton method failing when the previous point is not
in the immediate basin of attraction for the new point.
Presumably, this method has been developed by
John Hamal Hubbard
or John Milnor,
and it has been implemented by
Arnaud Chéritat,
Christian Henriksen
and
Mitsuhiro Shishikura
as well. It is used in Mandel since version 5.4, the current implementation
consists of QmnPlane::newtonRay() and
QmnPlane::rayNewton(). These were ported to
Java as well. The
paper
by Tomoki Kawahira
gives a detailed description and an error estimate for the approximation (1).

There are two approaches to draw the external ray of angle θ by computing
the external argument (3): Argument scanning means to check every pixel,
which is very slow, but might be improved by a recursive subdivision of the image.
Argument tracing means to construct a sequence of triangles along the ray.
(I have learned this technique from Chapter 14 in
this
book, where boundary tracing is used to draw escape lines.)
First two pixels at the border of the image are found, such that the
external argument is <θ at P_{0} and ≥θ at P_{1}.
Then compute the external argument at the third vertex P_{2}:

If it is <θ, then the ray is intersecting the edge from P_{1}
to P_{2}. Reflect the triangle along this edge to obtain the new
vertex P_{2}, rename the old vertex P_{2} to P_{0},
and keep P_{1}.

If it is ≥θ, then the ray is intersecting the edge from P_{0}
to P_{2}. Reflect the triangle along this edge to obtain the new
vertex P_{2}, rename the old vertex P_{2} to P_{1},
and keep P_{0}.
Repeat this process with the new triangle … and draw the ray along
the sequence of triangles. This basic idea can be improved in several ways:

To make the ray look smooth, draw line segments instead of single pixels.
Choose the way of reflection such that the distance to previously drawn
pixels is increased.

Check this distance also to avoid going in circles around the endpoint.
(This is not yet implemented in the current version of Mandel.)

When computing the external argument according to (3), adjust the principal
value of arg(z_{n}/(z_{n}c)) by ±2π, if
z_{n}=f_{c}^{n}(z) or
z_{n}=f_{c}^{n}(c) is in the triangle between 0, c,
and α_{c} as explained above.

Save the intermediate arguments arg(z_{n}/(z_{n}c)) for
the next pixel to check for jump discontinuities by comparison, and to
adjust them. (The previous technique shifts these closer to the boundary,
but does not remove them.) A drawback of this approach is that it restricts
the number of iterations to the length of the array, which means in
particular that the end of the ray may be lost.

When the image shows a very small part of the Julia set or Mandelbrot set,
maybe the starting point will not be found in spite of the previous
techniques. It can be searched for on the boundary of a virtual enlarged
image, but this will be impractical when it is necessary to enlarge the
image by several orders of magnitude.
Argument tracing has been used in Mandel since 1997, the current implementation
consists of mndlbrot::turnsign() and QmnPlane::traceRay().

Johannes Riedl has used a method, where the next point is chosen from a
finite collection of points at a prescribed distance from the last one.

The power series expansion of Φ_{c}^{1} or
Φ_{M}^{1} could be used as well to draw external
rays, but I expect it to converge extremely slowly close to the boundary.

Noting that the direction of an external ray is determined uniquely from the
fact that it is orthogonal to the equipotential line, external rays can be
drawn by solving a differential equation numerically. I am not sure if this
method will be numerically stable, and I have received different claims
about it.
The inverse of Φ_{c}(z) or Φ_{M}(c) can vary
considerably on a small interval: the angles 1/7 (of period 3) and
321685687669320/2251799813685247 (of period 51)
differ by about 10^{15}, but the landing points of the corresponding
parameter rays are about 0.035 apart. On this scale, the rays cannot be drawn
well by tracing the argument, unless highprecision
arithmetics is used. The following image shows the parameter ray for 1/7 and
two rays of period 51 in the slit betwen the main cardioid and the
1/3component. The argument is not distinguishable (left), but the
Newton method has no problem (right). In the dynamic plane, backwards iteration
works as well.
Note that most of the methods above can be modified to draw equipotential lines.
In this case there is no problem with the ambiguity of the argument. Boundary
scanning is slow and the Newton method has three problems: finding a starting
point, drawing several lines around a disconnected Julia set, and choosing
the number of points depending on the chosen subset of the plane. Therefore
I am using boundary tracing with starting points on several lines through the
image. See QmnPlane::green(). Still, sometimes not the whole line
is drawn, or not all components are found.
Last modified: January 2014.
Disclaimer: I am not responsible for linked sites by other people.