Rosetta Code

Mandelbrot set: Difference between revisions

Page Discussion

Mandelbrot set (view source)

Revision as of 08:23, 6 August 2022

657 bytes removed ,  1 year ago
→‎Normalized Counting, Distance Estimation, Mercator Maps and Perturbation Theory: Section final revised.
(→‎Normalized Counting, Distance Estimation, Mercator Maps and Perturbation Theory: Added reasons for choosing the rebasing condition abs(z) < abs(epsilon).)
Line 6,551:
savefig("Mandelbrot_set_3.png")</lang>
 
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.
<lang julia>using Plots
gr(aspect_ratio=:equal, axis=true, ticks=true, legend=false, dpi=200)
Line 6,578:
savefig("Mercator_Mandelbrot_map.png")
 
X, Y = real(C), imag(C) # zoom images (adjust circle size 6050 and zoom level 20 as needed)
R, c, z = 6050 .* 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)
Line 6,589:
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. 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. Placing the iteration in a separate function speeds up the calculation slightly.
<lang julia>using Plots
gr(aspect_ratio=:equal, axis=true, ticks=true, legend=false, dpi=200)
Line 6,607:
S = zeros(Complex{Float64}, n+1)
for k in 1:n+1
S[k] =if abs2(z) < abs2(r)
z = z ^ 2S[k] += cz
if abs2( z) >= z ^ 2 + abs2(r)c
else▼
error("n is too large (reference sequence diverges for reference point c)!")
end
Line 6,636 ⟶ 6,637:
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. In this example, moving the iteration to its own function results in a large speedup.
<lang julia>using Plots
gr(aspect_ratio=:equal, axis=true, ticks=true, legend=false, dpi=200)
Line 6,654 ⟶ 6,655:
S, dS = zeros(Complex{Float64}, n+1), zeros(Complex{Float64}, n+1)
for k in 1:n+1
S[k],if dS[k]abs2(z) = z,< dzabs2(r)
z, dz = z ^ 2 + cS[k], 2dS[k] *= z *, dz + 1
z, dz = z ^ 2 + c, 2 * z * dz + 1
if abs2(z) > abs2(r)▼
else
error("n is too large (reference sequence diverges for reference point c)!")
end
Line 6,670 ⟶ 6,672:
D = zeros(size(C))
 
iteration(S, dS, E, dE, C) = (2 .* S .+ E) .* E .+ C, 2 .* ((S .+ E) .* dE .+ EdS .* dSE)
for k in 1:n
M = abs2.(Z) .< abs2(r)
Line 6,683 ⟶ 6,685:
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.
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.▼
<lang julia>using Plots
gr(aspect_ratio=:equal, axis=true, ticks=true, legend=false, dpi=200)
Line 6,701 ⟶ 6,703:
S = zeros(Complex{Float64}, n+1)
for k in 1:n+1
S[k] =if abs2(z) < abs2(r)
z = z ^ 2S[k] += cz
if abs2( z) >= z ^ 2 + abs2(r)c
else
error("n is too large (reference sequence diverges for reference point c)!")
end
Line 6,714 ⟶ 6,717:
C = (.- 4.0) .* exp.((A' .+ B .* im) .* im)
 
Z, dZ, E = zero(C), zero(C), zero(C)
D = zeros(size(C))
 
iteration(s, epsilon, epsilondelta) = (2 * S[index]s + epsilon) * epsilon + delta▼
abs2_r = abs2(r)
for i in 1:h+1, j in 1:d+1
z, dz, epsilondelta = Z[i, j], dZ[i, j], EC[i, j]
deltaepsilon, index = C[i, j]z, 1
for k in 1:n
▲ if abs2(z) >< abs2(r)
▲ epsilon = (2 * S[index] + epsilon) * epsilon + delta
z, dz = S[index +if 1]abs2(z) +< abs2(epsilon,) # 2rebase *if z *is dzcloser +to 1zero
abs2_z epsilon, abs2_epsilonindex = abs2(z), abs2(epsilon)1
if abs2_z > abs2_r end
epsilon, index = ziteration(S[index], epsilon, 1delta)▼
indexz, dz = S[index + 1] + epsilon, 2 * z * dz + 1▼
index = index + 1
endelse▼
break
▲ end
if abs2_z < abs2_epsilon # rebasing if this reduces the distance to the reference sequence
▲ epsilon, index = z, 1
▲ else
▲ index = index + 1
end
end
Z[i, j], dZ[i, j], E[i, j] = z, dz, epsilon
end
 
Line 6,741 ⟶ 6,743:
 
heatmap(D' .^ 0.025, c=:nipy_spectral)
savefig("Mercator_Mandelbrot_base_mapMercator_Mandelbrot_rebase_map.png")</lang>
 
▲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, theThe 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) for 0 < gamma < 1 defines a circular area containing the origin. For gamma = 1/2 and below, the circle lies entirely within thea rescuebailout 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 haveare the same distanceequidistant from theO(0, origin0) and theS(6, reference point0). 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 theO(0, origin0) as from the referenceS(6, point0). All these points lie on a circle, withwhich has center M(-2, 0) and radius R = 4. This circle lies completely in a rescuebailout circle around the origin O(0, 0) with escape radius r = 6. If you choose a factor gamma > 1, you will get corresponding circles containing the reference point S. 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. 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
gr(aspect_ratio=:equal, axis=true, ticks=true, legend=false, dpi=200)
Line 6,776 ⟶ 6,778:
 
heatmap(D' .^ 0.025, c=:nipy_spectral)
savefig("Mercator_Mandelbrot_big_mapMercator_Mandelbrot_bigfloat_map.png")</lang>
 
=={{header|Kotlin}}==
Majow
305

edits

Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/105968"