Gesamtschule Aachen-Brand, 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
by Digia. Features include:
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
For your convenience, a Windows executable is offered in addition.
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 puzzle-pieces.
Display Julia sets on the Riemann
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.
by marking small Mandelbrot sets and small Julia sets.
embedded Julia sets
by a kind of preperiodic renormalization.
Display the asymptotic self-similarity at
on multiple scales, and the local similarity
to Julia sets.
Draw the parameter plane and Julia sets for other one-parameter
families, which satisfy critical relations or show a
persistent Siegel disk.
Available are unicritical polynomials, Branner-Fagella
polynomials, various families of cubic and quartic polynomials,
quadratic rational mappings,
transcendental functions, a simulation of quasiconformal surgery, the
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.
are providing an extensive introduction to basic and
advanced aspects of complex dynamics.
Start Mandel from the GAP package
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),
Download the source code (485 kB)
or the Windows exe (5.5 MB).
There is an executable for Mac as well
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:
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 DOS-program logistic.exe deals with the real iteration of z2 + c,
which is conjugate to the standard logistic map. It provides features like
drawing q-curves, finding hyperbolic intervals or periodic points, cobweb
graphics. This program may be used and redistributed without restriction.
It comes with no warranty.
2.0 of July 30, 1998. Preliminary manual logistic.ps.
How to draw external rays:
In the iteration of quadratic polynomials fc(z) = z2 + 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
The filled Julia set Kc contains those values z that are not
going to the
point at ∞
under iteration. The dynamics is chaotic on its boundary ∂Kc ,
which is the Julia set.
The mapping fc(z) = z2 + c is conjugate to
F(z) = z2 in a neighborhood of ∞. The conjugacy is provided
by the Boettcher mapping Φc(z), which satisfies
Φc(fc(z)) = F(Φc(z)). We have
(1) Φc(z) = limn→∞
for a suitable choice of the 2n-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/(fcn-1(z))2)
= z ∏n=1∞
In the dynamics of quadratic polynomials, there is a basic dichotomy:
The mapping ΦM(c) is defined by evaluating the Boettcher mapping
at the critcal value z = c, i.e., ΦM(c) = Φc(c).
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
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 well-defined at the
critcal value z = c.
Computation of the external argument:
External rays are characterized by the fact that the argument argc(z) of
Φc(z), or the argument argM(c) of ΦM(c),
equals a given value θ. For large |z|, formula (2) implies
(3) argc(z) = arg(z) + ∑n=1∞
1/2n arg(1 + c/(fcn-1(z))2)
= arg(z) + ∑n=1∞1/2n
Here arg(z) denotes the
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 single-valued. 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 zn=fcn(z) of
z and the arguments arg(zn/(zn-c))
are computed. The function arg(z/(z-c)) 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/(z-c) 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:
When the modified function arg(z/(z-c)) is used instead of the principal value
for computing argc(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 zn is in
the triangle between 0, c, and αc, adjust the principal value
of arg(zn/(zn-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 argc(z), and thus, external rays.
The computed value of argc(z), and thus the color, is jumping on
certain curves joining preimages of 0, where some iterate zn is
crossing the line where arg(z/(z-c)) is discontinuous. In the left image, the
principal value of arg(z/(z-c)) is used, and in the right image, the argument
is adjusted. The erroneous discontinuities of the computed argc(z)
are moved closer to the Julia set.
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/(z-c))
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/(z-c)) 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.
The same modification seems to work for the computation of the argument
argM(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/(z-c)) 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
See also the algorithms 7 and 8, where no discontinuity appears, but only the
ends of certain rays are visible.
has suggested that boundary scanning (below) can be used to draw several rays simultaneously.
Drawing a single ray:
The following methods apply with small modifications to the dynamic plane and
the parameter plane as well:
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 zl,j corresponds to a radius
R1/2l and the angle 2jθ. Then
fc(z) maps zl,j to z(l-1),(j+1). This point,
which was constructed before, has two preimages under fc(z) . The
one that is closer to z(l-1),j is the correct one.
This criterion was proved by
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 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 high-precision
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/3-component. The argument is not distinguishable (left), but the
Newton method has no problem (right). In the dynamic plane, backwards iteration
works as well.
Choose a sequence of radii tending to 1, e.g.,
Rn = R1/2n with R large.
Search corresponding points on the ray, which satisfy
Φc(zn) = Rn eiθ or
ΦM(cn) = Rn eiθ,
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) = fcn(z) - R ei2nθ,
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
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
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
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 P0 and ≥θ at P1.
Then compute the external argument at the third vertex P2:
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:
If it is <θ, then the ray is intersecting the edge from P1
to P2. Reflect the triangle along this edge to obtain the new
vertex P2, rename the old vertex P2 to P0,
and keep P1.
If it is ≥θ, then the ray is intersecting the edge from P0
to P2. Reflect the triangle along this edge to obtain the new
vertex P2, rename the old vertex P2 to P1,
and keep P0.
Argument tracing has been used in Mandel since 1997, the current implementation
consists of mndlbrot::turnsign() and QmnPlane::traceRay().
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(zn/(zn-c)) by ±2π, if
zn=fcn(c) is in the triangle between 0, c,
and αc as explained above.
Save the intermediate arguments arg(zn/(zn-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.
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
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.