Mandelbrot set: Difference between revisions
Content deleted Content added
→Normalized Counting, Distance Estimation, Mercator Maps and Perturbation Theory: Added reasons for choosing the rebasing condition abs(z) < abs(epsilon). |
→Normalized Counting, Distance Estimation, Mercator Maps and Perturbation Theory: Section final revised. |
||
Line 6,551: | Line 6,551: | ||
savefig("Mandelbrot_set_3.png")</lang> |
savefig("Mandelbrot_set_3.png")</lang> |
||
A small change in the code above creates Mercator maps and zoom images of the Mandelbrot set. The |
A small change in the code above creates Mercator maps and zoom images of the Mandelbrot set. The abs2 function is used to speed up the calculations inside the loop. See also the album [https://www.flickr.com/photos/arenamontanus/albums/72157615740829949 Mercator Mandelbrot Maps] by Anders Sandberg. |
||
<lang julia>using Plots |
<lang julia>using Plots |
||
gr(aspect_ratio=:equal, axis=true, ticks=true, legend=false, dpi=200) |
gr(aspect_ratio=:equal, axis=true, ticks=true, legend=false, dpi=200) |
||
Line 6,578: | Line 6,578: | ||
savefig("Mercator_Mandelbrot_map.png") |
savefig("Mercator_Mandelbrot_map.png") |
||
X, Y = real(C), imag(C) # zoom images (adjust circle size |
X, Y = real(C), imag(C) # zoom images (adjust circle size 50 and zoom level 20 as needed) |
||
R, c, z = |
R, c, z = 50 .* 2 ./ d .* pi .* ones(d+1)' .* exp.(.- B), min(d, h) + 1, max(0, h - d) ÷ 20 |
||
gr(aspect_ratio=:equal, axis=false, ticks=false, legend=false, markerstrokewidth=0, dpi=200) |
gr(aspect_ratio=:equal, axis=false, ticks=false, legend=false, markerstrokewidth=0, dpi=200) |
||
Line 6,589: | Line 6,589: | ||
savefig("Mercator_Mandelbrot_plot.png")</lang> |
savefig("Mercator_Mandelbrot_plot.png")</lang> |
||
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 |
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. Placing the iteration in a separate function speeds up the calculation. 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. |
||
<lang julia>using Plots |
<lang julia>using Plots |
||
gr(aspect_ratio=:equal, axis=true, ticks=true, legend=false, dpi=200) |
gr(aspect_ratio=:equal, axis=true, ticks=true, legend=false, dpi=200) |
||
Line 6,607: | Line 6,607: | ||
S = zeros(Complex{Float64}, n+1) |
S = zeros(Complex{Float64}, n+1) |
||
for k in 1:n+1 |
for k in 1:n+1 |
||
if abs2(z) < abs2(r) |
|||
S[k] = z |
|||
z = z ^ 2 + c |
|||
⚫ | |||
error("n is too large (reference sequence diverges for reference point c)!") |
error("n is too large (reference sequence diverges for reference point c)!") |
||
end |
end |
||
Line 6,636: | Line 6,637: | ||
savefig("Mandelbrot_deep_zoom.png")</lang> |
savefig("Mandelbrot_deep_zoom.png")</lang> |
||
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 |
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. |
||
<lang julia>using Plots |
<lang julia>using Plots |
||
gr(aspect_ratio=:equal, axis=true, ticks=true, legend=false, dpi=200) |
gr(aspect_ratio=:equal, axis=true, ticks=true, legend=false, dpi=200) |
||
Line 6,654: | Line 6,655: | ||
S, dS = zeros(Complex{Float64}, n+1), zeros(Complex{Float64}, n+1) |
S, dS = zeros(Complex{Float64}, n+1), zeros(Complex{Float64}, n+1) |
||
for k in 1:n+1 |
for k in 1:n+1 |
||
if abs2(z) < abs2(r) |
|||
S[k], dS[k] = z, dz |
|||
z, dz = z ^ 2 + c, 2 * z * dz + 1 |
|||
⚫ | |||
else |
|||
error("n is too large (reference sequence diverges for reference point c)!") |
error("n is too large (reference sequence diverges for reference point c)!") |
||
end |
end |
||
Line 6,670: | Line 6,672: | ||
D = zeros(size(C)) |
D = zeros(size(C)) |
||
iteration(S, dS, E, dE, C) = (2 .* S .+ E) .* E .+ C, 2 .* ((S .+ E) .* dE .+ |
iteration(S, dS, E, dE, C) = (2 .* S .+ E) .* E .+ C, 2 .* ((S .+ E) .* dE .+ dS .* E) |
||
for k in 1:n |
for k in 1:n |
||
M = abs2.(Z) .< abs2(r) |
M = abs2.(Z) .< abs2(r) |
||
Line 6,683: | Line 6,685: | ||
savefig("Mercator_Mandelbrot_deep_map.png")</lang> |
savefig("Mercator_Mandelbrot_deep_map.png")</lang> |
||
Another approach to reducing glitches is 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. |
|||
⚫ | |||
<lang julia>using Plots |
<lang julia>using Plots |
||
gr(aspect_ratio=:equal, axis=true, ticks=true, legend=false, dpi=200) |
gr(aspect_ratio=:equal, axis=true, ticks=true, legend=false, dpi=200) |
||
Line 6,701: | Line 6,703: | ||
S = zeros(Complex{Float64}, n+1) |
S = zeros(Complex{Float64}, n+1) |
||
for k in 1:n+1 |
for k in 1:n+1 |
||
if abs2(z) < abs2(r) |
|||
S[k] = z |
|||
z = z ^ 2 + c |
|||
else |
|||
error("n is too large (reference sequence diverges for reference point c)!") |
error("n is too large (reference sequence diverges for reference point c)!") |
||
end |
end |
||
Line 6,714: | Line 6,717: | ||
C = (.- 4.0) .* exp.((A' .+ B .* im) .* im) |
C = (.- 4.0) .* exp.((A' .+ B .* im) .* im) |
||
Z, dZ |
Z, dZ = zero(C), zero(C) |
||
D = zeros(size(C)) |
D = zeros(size(C)) |
||
⚫ | |||
abs2_r = abs2(r) |
|||
for i in 1:h+1, j in 1:d+1 |
for i in 1:h+1, j in 1:d+1 |
||
z, dz, |
z, dz, delta = Z[i, j], dZ[i, j], C[i, j] |
||
epsilon, index = z, 1 |
|||
for k in 1:n |
for k in 1:n |
||
⚫ | |||
⚫ | |||
if abs2(z) < abs2(epsilon) # rebase if z is closer to zero |
|||
epsilon, index = z, 1 |
|||
end |
|||
⚫ | |||
⚫ | |||
index = index + 1 |
|||
⚫ | |||
break |
break |
||
⚫ | |||
if abs2_z < abs2_epsilon # rebasing if this reduces the distance to the reference sequence |
|||
⚫ | |||
⚫ | |||
⚫ | |||
end |
end |
||
end |
end |
||
Z[i, j], dZ |
Z[i, j], dZ[i, j] = z, dz |
||
end |
end |
||
Line 6,741: | Line 6,743: | ||
heatmap(D' .^ 0.025, c=:nipy_spectral) |
heatmap(D' .^ 0.025, c=:nipy_spectral) |
||
savefig(" |
savefig("Mercator_Mandelbrot_rebase_map.png")</lang> |
||
⚫ | The condition abs(z) < abs(epsilon) does not define a circular area but a half-plane. If you insert a factor, the condition abs(z) < gamma * abs(epsilon) for 0 < gamma < 1 defines a circular area containing the origin. For gamma = 1/2 and below, the circle lies entirely within a bailout circle around the origin with escape radius r > 2. For example, if the reference point is at S(6, 0), the equation abs(z) = abs(epsilon) defines the set of all points that are equidistant from O(0, 0) and S(6, 0). 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 O(0, 0) as from S(6, 0). All these points lie on a circle, which has center M(-2, 0) and radius R = 4. This circle lies completely in a bailout circle around the origin O(0, 0) with escape radius r = 6. If you choose a factor gamma > 1, you get corresponding circles containing the reference point S. However, the condition abs(z) < abs(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. You can only rebase to the beginning of the reference sequence. If you want to rebase to the reference point S[m], the distance epsilon = z - S[m] must be calculated. With this subtraction, however, the entire precision is lost due to cancellation. Only at the first reference point there is no cancellation due to S[1] = 0, because epsilon = z - S[1] = z - 0 = z. The images can be verified by a (slow) calculation with BigFloats. |
||
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. |
|||
<lang julia>using Plots |
<lang julia>using Plots |
||
gr(aspect_ratio=:equal, axis=true, ticks=true, legend=false, dpi=200) |
gr(aspect_ratio=:equal, axis=true, ticks=true, legend=false, dpi=200) |
||
Line 6,776: | Line 6,778: | ||
heatmap(D' .^ 0.025, c=:nipy_spectral) |
heatmap(D' .^ 0.025, c=:nipy_spectral) |
||
savefig(" |
savefig("Mercator_Mandelbrot_bigfloat_map.png")</lang> |
||
=={{header|Kotlin}}== |
=={{header|Kotlin}}== |