Mandelbrot set: Difference between revisions

A small change in the code above creates Mercator maps and zoom images of the Mandelbrot set. The squared abs2 function is used to speed up the calculations inside the loop because it avoids a time-consuming square root calculation. However, the derivative dZ can be very large, so that an overflow error can occur when calculating abs2(dZ). Therefore, when calculating the distance, the abs function is used again, which avoids such errors due to a clever implementation of the hypot function. See also the album [https://www.flickr.com/photos/arenamontanus/albums/72157615740829949 Mercator Mandelbrot Maps] by Anders Sandberg.

X, Y = real(C), imag(C) # zoom images (adjust circle size 60 and zoom level 20 as needed)

R, c, z = 60 .* 2 ./ d .* pi .* ones(d+1)' .* exp.(.- B), min(d, h) + 1, max(0, h - d) ÷ 20

For deep zoom images it is sufficient to calculate a single point with high accuracy. A good approximation can then be found for all other points by means of a perturbation calculation with standard accuracy. See [https://en.wikipedia.org/wiki/Plotting_algorithms_for_the_Mandelbrot_set#Perturbation_theory_and_series_approximation Perturbation theory] (Wikipedia) and [https://gbillotey.github.io/Fractalshades-doc/math.html Mathematical background] (Fractalshades) for more details. Placing the iteration in a separate function speeds up the calculation slightly.

S[k] = z

z = z ^ 2 + c

if abs2(z) > abs2(r)

Of course, deep Mercator maps can also be created. See also the image [https://www.flickr.com/photos/arenamontanus/3430921497/in/album-72157615740829949/ Deeper Mercator Mandelbrot] by Anders Sandberg. To reduce glitches, an additional reference sequence for the derivation (dS) is recorded with high precision. In this example, moving the iteration to its own function results in a large speedup.

S[k], dS[k] = z, dz

z, dz = z ^ 2 + c, 2 * z * dz + 1

if abs2(z) > abs2(r)

iteration(S, dS, E, dE, C) = (2 .* S .+ E) .* E .+ C, 2 .* ((S .+ E) .* dE .+ E .* dS)

Another approach to reduce the glitches is the so-called ''rebasing''. See [https://gbillotey.github.io/Fractalshades-doc/math.html#avoiding-loss-of-precision Avoiding loss of precision] (Fractalshades) and [https://fractalforums.org/fractal-mathematics-and-new-theories/28/another-solution-to-perturbation-glitches/4360 Another solution to perturbation glitches] (Fractalforums) for details. Remarkably, the condition for the rebasing abs2(z) < abs2(epsilon) or abs(z) < abs(epsilon) does not define a circular area but a half-plane. If you want a circular area, you can insert a factor. For 0 < gamma < 1 the condition abs(z) < gamma * abs(epsilon) defines a circular area containing the origin. For gamma = 1/2 and below, the circle lies entirely within the rescue circle around the origin O(0, 0) with escape radius r > 2, even if the associated reference point S is on the edge of the rescue circle. For example, if the reference point is at S(6, 0), the equation abs(z) = abs(epsilon) defines the set of all points that have the same distance from the origin and the reference point. This is exactly the perpendicular bisector x = 3. Considering instead the equation abs(z) = 1/2 * abs(epsilon), one finds all points half as far from the origin as from the reference point. All these points lie on a circle with center M(-2, 0) and radius R = 4. This circle lies completely in a rescue circle around the origin with escape radius r = 6. If you choose a factor gamma > 1, you will get corresponding circles containing the reference point. However, the condition abs(z) < abs(epsilon) or abs2(z) < abs2(epsilon) is optimal because it keeps the differences epsilon to the reference sequence as small as possible. Rebasing occurs if and only if the distance abs(epsilon) to the reference sequence becomes smaller as a result.

S[k] = z

z = z ^ 2 + c

if abs2(z) > abs2(r)

Z, dZ, E = zero(C), zero(C), zero(C)

abs2_r = abs2(r)

z, dz, epsilon = Z[i, j], dZ[i, j], E[i, j]

delta, index = C[i, j], 1

epsilon = (2 * S[index] + epsilon) * epsilon + delta

z, dz = S[index + 1] + epsilon, 2 * z * dz + 1

abs2_z, abs2_epsilon = abs2(z), abs2(epsilon)

if abs2_z > abs2_r

end

if abs2_z < abs2_epsilon # rebasing if this reduces the distance to the reference sequence

epsilon, index = z, 1

else

index = index + 1

Z[i, j], dZ[i, j], E[i, j] = z, dz, epsilon

savefig("Mercator_Mandelbrot_base_map.png")</lang>

The images can be verified by a (slow) calculation with BigFloats. There are libraries that are faster than BigFloats, for example DoubleFloats.jl and MultiFloats.jl. Unfortunately, the precision of the DoubleFloats is not sufficient, and only the double MultiFloats (Float64x2) are fast. For deep zoom images you need at least triple or quadruple MultiFloats (Float64x3, Float64x4), but these are slower than BigFloats.

savefig("Mercator_Mandelbrot_big_map.png")</lang>