Wolf Jung   jung@mndynamics.com
Gesamtschule Aachen-Brand, Rombachstrasse 41, 52078 Aachen, Germany.

Mandel: software for real and complex dynamics

• Draw the Mandelbrot set and Julia sets in two windows at the same time, in one window, or in full-screen mode. Available methods of drawing include escape time, hyperbolic components or attracting basins, distance estimate, marking external rays or puzzle-pieces.
• 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.
• Find centers and Misiurewicz points in the parameter plane, or periodic and preperiodic points in the dynamic plane. Illustrate renormalization by marking little Mandelbrot sets and little Julia sets.
• Display the asymptotic self-similarity at Misiurewicz points 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, Newton's method, transcendental functions, and a simulation of quasi-conformal surgery.
• Save images as PostScript *.eps, save or load a *.png image.
• Future plans include an extensive demo and manual, the spider algorithm, the construction of Hubbard trees, and incorporating real iteration.

The non-GUI-version has been renamed: Exmandel is a DOS-program with VGA-graphics. Its mathematical functions are similar to Mandel, but it is controlled by keyboard only. This program may be used and redistributed without restriction. It comes with no warranty.
Download exmandel.exe 4.1 of June 21, 2007 (208 kB).

Other sites with programs on complex dynamics or fractals:
fractint, IMS, Arnaud Chéritat, Xaos, spanky.

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 z² + 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.
Download logistic.exe 2.0 of July 30, 1998. Preliminary manual logistic.ps.

How to draw external rays:

Background:

(1)  Φ_c(z) = lim_n→∞ (f_cⁿ(z))^1/2n

for a suitable choice of the 2ⁿ-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^n-1(z))²) ^1/2n = z ∏_n=1^∞ (f_cⁿ(z)/(f_cⁿ(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 well-defined at the critcal value z = c.

Computation of the external argument:

(3)  arg_c(z) = arg(z) + ∑_n=1^∞ 1/2ⁿ arg(1 + c/(f_c^n-1(z))²) = arg(z) + ∑_n=1^∞1/2ⁿ arg(f_cⁿ(z)/(f_cⁿ(z)-c)).

Note that the principal value of the argument satisfies -π<arg(z)≤π, and it is discontinuous on the negative real axis. The function arg(z/(z-c)) is discontinuous on the straight line segment from z = 0 to z = c. When the Julia set is connected, the discontinuity can be moved such that the argument is continuous in the whole exterior, and (3) is valid everywhere. Approximating the Julia set by two straight line segments from the fixed point α_c to 0 and c, respectively, the argument can be adjusted by three simple inequalities, cf. mndlbrot::turn(). The same formulas work in the parameter plane as well. The left image shows the discontinuities arising from the principal value, 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:

Drawing a single ray:

• In the dynamic plane, external rays can be drawn by backwards iteration: a point z with argument θ is a preimage of a point with argument 2θ under f_c(z). The criterion for choosing the correct preimage was proved by Thierry Bousch. The following methods apply with small modifications to the dynamic plane and the parameter plane as well:
• There are two approaches to draw the external ray of angle θ by computing the external argument (3): Boundary scanning means to check every pixel, which is very slow, but might be improved by a recursive subdivision of the image. Boundary tracing means to construct a sequence of triangles along the ray. (I have learned this technique from Chapter 14 in this book, where it 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₀ and ≥θ at P₁. Then compute the external argument at the third vertex P₂:
  • If it is <θ, then the ray intersects the edge from P₁ to P₂. Reflect the triangle along this edge to obtain the new vertex P₂, rename the old vertex P₂ to P₀, and keep P₁.
  • If it is ≥θ, then the ray intersects the edge from P₀ to P₂. Reflect the triangle along this edge to obtain the new vertex P₂, rename the old vertex P₂ to P₁, and keep P₀.
  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, adjust the principal value if f_cⁿ(z) is in the triangle between 0, c, and α_c as explained above.
  • Save the intermediate arguments from the previous 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.
  Boundary tracing has been used in Mandel since 1997, the current implementation consists of mndlbrot::turnsign() and QmnPlane::exray(). I shall soon implement the following method 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, Newton's 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ⁿ(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 Newton's 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 and Christian Henriksen as well.
• 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 extremly 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 from solving a differential equation numerically. I am not sure if this method will be numerically stable.