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 The Qt Company. 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 provide an extensive introduction to basic and advanced aspects of complex dynamics.
 Start Mandel from the GAP package IMG 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).
Download the source code (0.5 MB) or the Windows exe (5.7 MB).
Note: For polynomial spiders and twists, you need an old version in addition.
There is an executable for Mac as well (7.6 MB, Mandel 5.10, no matings).
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, Hiroyuki Inou, Chris King, 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.
Download logistic.exe 2.0 of July 30, 1998. Preliminary manual logistic.ps.
How to draw external rays:
Show more ...
Hide again ...
(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:
(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:
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.
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.
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.
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.
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().
 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}.
 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.
 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.
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.