Julia sets: True shape and escape time

  • 57 Replies
  • 1752 Views

0 Members and 1 Guest are viewing this topic.

Offline marcm200

  • *
  • Fractal Flamingo
  • ****
  • Posts: 340
« Reply #15 on: May 21, 2019, 10:41:51 AM »
Finding a tight bounding box for the sectionally defined Julia set is tricky.

I worked on implementing a simple rule for the quadratic case: if imag(z)<0 use c1 otherwise use c2.

But that separation type got me thinking. The x-axis is always touched by two adjacent pixels and their underlying squares, one where the x axis is the top edge, one where it is the bottom edge. So in one of those two pixels the bounding box for f(horizontal edges) will be computed with two different c values and the usual bounding box that calculates the images of one iteration of the corners of the pixel is invalid.

I settled with a very conservative (basically un-informative) approach: If a square touches the x-axis, its bounding box under f is constantly set to the whole [-2..+2] x [-2..+2] square. That's certainly correct - but does not allow the propagation of any color.

I'm thinking of something like this is a better and informative way: calculating the bounding boxes for such a problematic square with the c value fixed and use both. And then use a rectangle that encompasses both individual bounding boxes. But I cannot give a mathematical reasoning for that. Any ideas?

What I cannot explain in the image below is the occurances of quite a lot of white islands in gray areas that are surrounded by lots of black (upper left quadrant, marked with a red circle). So I wonder whether in higher resolutions those will eventually all form one big white area - or if the islands persists.

The section definition imag(z)<0 introduces a non-continuous jump into the function. Is that irregularity causing the issue? Does anyone have an explanation for that type of behaviour?

And those protruding sticks/arcs look a bit like mappings from the x-axis to other parts. I hope they can be educed and set to a color if I figured out an informative bounding box.

Question for the moderators: Since my results are mostly images (and some polygon point lists): Should the thread be moved to the image thread section of the forum?




Offline claude

  • *
  • 3e
  • *****
  • Posts: 1149
    • mathr.co.uk
« Reply #16 on: May 21, 2019, 11:25:06 AM »
Walk the plane, compute individual Julia sets for the quadratic case (hence any Mandelbrot image pixel corresponds to one single point in the complex plane, not a square like for the Julia sets) - only the white parts and a small region around the origin were actually computed - and check whether the origin can be colored black or white.

You should use an interval for Mandelbrot C coordinates when computing the Julia sets, otherwise you can't trust the image.

Offline marcm200

  • *
  • Fractal Flamingo
  • ****
  • Posts: 340
« Reply #17 on: May 21, 2019, 11:54:57 AM »
You should use an interval for Mandelbrot C coordinates when computing the Julia sets, otherwise you can't trust the image.
Yes, that's the trustworthy way for a complete view of the Mandelbrot set. But my final goal was (and still is) to get a list of complex numbers with rational components up to a certain digit length. And for that purpose it was sufficient to use point sampling - which proved very very slow anyways. I wouldn't want to imagine how long a truely trustworthy Mandelbrot image would take (it's definitely faster with the reverse cell graph, but as the authors stated: That's way more complicated.)



Offline claude

  • *
  • 3e
  • *****
  • Posts: 1149
    • mathr.co.uk
« Reply #18 on: May 21, 2019, 12:59:36 PM »
I'd be interested to see what happens with the Burning Ship, both the parameter plane and its (quasi-)Julia sets.  The Burning Ship has some areas that are neither escaping to infinity nor attracted to a periodic cycle - curious to see how they will be coloured.  Some interesting example C for Julia sets include (-3/4,-1/2), (1/4,-3/4)

Offline marcm200

  • *
  • Fractal Flamingo
  • ****
  • Posts: 340
« Reply #19 on: May 22, 2019, 08:25:21 AM »
As claude suggested, I tried the Burning Ship formula for the two Julia set C values he posted (the first two in the image below).

Those with interior were quickly generated (only minutes in a spatial resolution of 2-10) with black found pretty late but than in an ever faster manner.

@claude: Are there specific regions in those Julia sets that show that non-period-but-bounded behaviour? (The origin?) I'm implementing a tool to show - for a starting pixel - the set of pixels it visits until it's stable, so one could have a look on where those orbits jump around in the plane.

As for the parameter plane, currently I cannot compute trustworthy images. I started reading about interval arithmetics and plan to code that.

The bounding box of f(A) is the same as for then standard Julia set. Here's my line of reasoning:

Let zn=x+i*y. x,y, real numbers, c=e+d*i, e,d real numbers.

Then zn+1 = (|x|+i*|y|)2 + (e+d*i) = (x▓-y▓+e) + i*(2*|x*y|+d)

For x,y both positive or both negative (or zero, since it does not discontinuously jump as the imag(z) rule of one of the sectionally defined sets), one can omit the abs function and the formula is the same as for the standard quadratic Julia set case.

For x negative and y positive (or the other way around), 2i|xy|+d = i*(-2*x*y+d), hence bracket part is a normal polynomial for which the same argument as in the paper can be used: it's a monotonic function since the grid of the plane is spaced in a way that no pixel  (x0..x1,y0..y1) ever crosses one of the coordinate axis. It touches but in a continuous way.

So the bounding box f(A) is the same as for the standard quadratic case.

Although by no means a proof, I did a computational check (16k x 16k grid of the -2..+2 square, each grid point was thought as a square A, I calculated the bounding box and checked 256 equally spaced points within A to see whether their image under f is in the computed bounding box. 68 trillion points were tested and passed. So the bounding box is not completely wrong. (Everything is a fraction of a power of two with an enumerator being a whole number so exact calculation for one iteration was possible with the usual C++ double precision data type).





Offline claude

  • *
  • 3e
  • *****
  • Posts: 1149
    • mathr.co.uk
« Reply #20 on: May 22, 2019, 05:28:57 PM »
Thanks!

@claude: Are there specific regions in those Julia sets that show that non-period-but-bounded behaviour? (The origin?) I'm implementing a tool to show - for a starting pixel - the set of pixels it visits until it's stable, so one could have a look on where those orbits jump around in the plane.
Yes the origin seems to go to an aperiodic bounded orbit for the C values I mentioned.  See https://math.stackexchange.com/questions/3234950/bounded-aperiodic-orbits-in-the-burning-ship  I think the trustworthy algorithm can provide part of the answer ("it remains bounded") though proving it is aperiodic might be more difficult?

Offline marcm200

  • *
  • Fractal Flamingo
  • ****
  • Posts: 340
« Reply #21 on: May 24, 2019, 09:47:43 AM »
I followed the orbit of the origin in a trustworthy manner - for claude's C values/burning ship and some standard quadratic ones.

I calculated in which pixel 0+0i lies at the lower left corner, and conservatively took this and its 8 neighbours as a starting set. The bounding box under one iteration of each of those pixels was conservatively calculated and then the ones from those pixels and so on until no unvisited pixels were hit.

The red regions in the images are those pixels that keep to themselves under f, so they are invariant.

The first image contains the burning ship Julia sets, all non-trustworthily downsized by irfanview, but the overall features stayed intact. Standard spatial resolution was 2-10 (4096 pixel, range -2..+2 in both directions). Origin is in the center of the image.

The second row shows an interesting feature. The big red region does not contain the origin. However if one looks at the 4096 version there are tiny red pixels in the image center. I attached the bigger image as a zip file.

And as comparison I calculated a couple of standard quadratic Julia sets.

If someone knows of interesting seed values of connected sets, I would be interested in checking them - provided they will produce black pixels (e.g. high periodicity or parabolic ones, aperiodic ones here as well?).

Is the hole in the first burning ship set a clue for an aperiodic orbit? I haven't seen that behaviour in the standard ones.

Or is the way the red region changes (or does not) when increasing the spatial resolution a way for an answer? The pixel represent ever smaller getting squares - but I don't now "how much" of a square in the bounding box is actually hit by an orbit, so I can't tell if it is safe to assume the red regions get smaller if one approaches a point cycle.



Offline marcm200

  • *
  • Fractal Flamingo
  • ****
  • Posts: 340
« Reply #22 on: May 26, 2019, 01:41:10 PM »
I had an idea about the black search algorithm to not have to follow all paths from one square to see whether they lie in a set of gray pixels that is invariant under f. If - in the way of following the paths for a starting pixel A a pixel B was hit, that is gray and was visited and followed in an earlier round of black-checking but after the last change to the image, then B was not part of an invariant set, so A cannot be either for now and the construction process of the visitable pixels for A can be stopped with the result "not able to be colored black for now". Depending on "how much" of the set still is out there, that could save a lot of time.

The good thing here is, this (maybe premature) stop cannot produce falsely marked black pixels. So even if the idea is wrong or the implementation faulty, in a sense it does not matter because even if I judge a pixel as "unknown" that might be colored black if the set was built completely, the method at most produces the same pixels as without the new routine. So the images remain trustworthy: black means bounded orbits, white means escaping ones and gray contains the boundary however disconnected it is.

(This idea is probably described in the paper or an intrinsic property of the reverse cell graph (which I'm not using), but then I didn't realize it).

That sped things up tremendously for some images (for c=-1 in the quadratic case about 200-fold). I checked all the parameters in the paper and luckily found that the optimized algorithm produced the same outcomes as before (I still cannot find black for the parabolic case c=0.25, the managable resolution of 2^-8 is too low).

Now it was possible to compute the image in the burning ship post (c=-0.625+0.6875i, standard quadratic case) to 4-times the spatial resolution (2-12) in about half an hour completely (I didn't like having to stop the computation early). And it showed more black and smaller red regions for the path of the origin that seem to become disconnected.

I think I need to buy more RAM now (the paper says the limit in 2013 was about 2-20 - quite some way to go there.










Offline marcm200

  • *
  • Fractal Flamingo
  • ****
  • Posts: 340
« Reply #23 on: May 28, 2019, 12:43:54 PM »
I incorporated some higher degree polynomials into the TSA.

After spending a day trying to figure out the bounding box for z3+a*z+c by hand using intervall arithmetics and always getting the subscripts of my intervals [x0..x1] and [y0..y1] (z=x+i*y) mixed up, I decided it's faster and easier to write a program that does symbolic interval arithmetics for the operations I need.

So adding new polynomials now is quite simple (I mainly do copying&pasting): WolframAlpha calculates the component functions for the real and imaginary zn+1, I put those into my IA parser that calculates the bounding box to a very un-optimized way (lots of max and min functions, lots of repeating powers at different positions etc.) and put that into the TSA (optimization did not lead to speed gains, I guess today's compilers are so efficient, they can do more than I can.)

Here are some examples (axis range was adjusted following the formula by Douady and is usually -3..+3 or greater). The bounding boxes I use are in the attached text file (zipped, too long to display here).

I've put one "Julia set" there that is all but trustworthy. It's one I calculated when I manually computed the bounding box. But the shape lead immeadiately to the thought: "That's not right" (which it wasn't).

CAUTION: The bounding boxes for order 6 and 8 polynomials contain terms like x8-y8 or x6+y6 - and depending on where in the plane x and y lie (x might be 2-12, y might be 20), one might/will run into rounding errors. So those images are more an example of how complicated the bounding boxes can get but still they can be evaluated pretty fast (arbitrary precision or integer arithmetics would fix that, but that's for the farther future on my list).

I wonder whether the TSA could also handle rational or transcendental functions (I think bounding boxes are universally working). There is no mentioning of that in the paper.

But now I'll focus on constructing polygons for the black region(s).




Offline claude

  • *
  • 3e
  • *****
  • Posts: 1149
    • mathr.co.uk
« Reply #24 on: May 28, 2019, 05:15:00 PM »
+1 for symbolic evaluation, bonus points if you can generate source code (or even binary code in memory that you can execute right away) directly from the input formula (I do it for perturbation and derivatives in my 'et' project)

Rational/transcendental functions may end up without any white regions, only black and grey, because escape methods are not appropriate for them (arbitrarily large values can become small again, depending on the function).  I haven't checked if the algorithm still works without white, maybe there will be an issue about finding a bounding box for the whole Julia set because rendering just a part won't work?  There is also the problem of multiple critical points if you want to render the parameter space and not just Julia sets.

Offline marcm200

  • *
  • Fractal Flamingo
  • ****
  • Posts: 340
« Reply #25 on: May 29, 2019, 07:13:11 AM »
+1 for symbolic evaluation, bonus points if you can generate source code (or even binary code in memory that you can execute right away) directly from the input formula (I do it for perturbation and derivatives in my 'et' project)
Thanks. Directly generated binary code - neat! But that's beyond my skill set. The last time I worked with binary/assembler was in the good old Amstrad days. My parser outputs just a text file whose content I can copy&paste into the source from my implementation of the TSA, but I have to compile it every time I add a new bounding box. So no "user-defined formula" feature here.

Rational/transcendental functions may end up without any white regions, only black and grey, because escape methods are not appropriate for them (arbitrarily large values can become small again, depending on the function).  I haven't checked if the algorithm still works without white, maybe there will be an issue about finding a bounding box for the whole Julia set because rendering just a part won't work? 
That's a fair point. If nothing escapes, then the trilogy of bounded, escaping and the gray area containing the boundary does not apply. Maybe the method could detect some sort of functional cycle of the sort: sin(sin(sin(square A)...) is a proper subset of [-1..1].

The paper says their algorithm is global in nature and needs always a complete region that definitely cirumferences the Julia set of interest. So zooming is not possible. I think it's because of the rules "all paths lead to" and "no path leads to", where one inherently needs to know "everything", or at least everything that can be visited from one starting square A. So if one would only need that color, it might be possible to work around checking the complete plane (I put that on my to-do-list as well).
 

Offline pauldelbrot

  • *
  • 3f
  • ******
  • Posts: 1292
« Reply #26 on: May 29, 2019, 03:28:41 PM »
For rational maps, I'd suggest the following procedure would work:

  • Follow all critical orbits for some large number of iterations.
  • For each one: Run it a few thousand more iterations and see if it gets close to the landing points of any previous ones. If it does, discard it.
  • Divide the Riemann sphere into small regions in a grid; e.g. by keeping two square bitmaps, representing in one the inside of the unit circle and in the other the outside of it by representing the inside of the unit circle for w = 1/z.
  • Color the pixels in these bitmaps containing the remaining critical orbit landing points red, green, blue, etc.; maybe use a conservative distance estimate at each one to widen this to a small disk of pixels.
  • Apply the original algorithm, only instead of looking for pixels to hit white or not, you look for them to hit all those various colors, and propagate those colors.

In the end you should have colorful filled-in basins of attraction along with grey along the Julia set itself. Conversion of the two bitmaps into a suitable visual representation of a sphere is left as an exercise for the reader.

Offline marcm200

  • *
  • Fractal Flamingo
  • ****
  • Posts: 340
« Reply #27 on: May 30, 2019, 01:45:47 PM »
@pauldelbrot: And here I was, thinking that adding a function and a bounding box would suffice... Thanks for the detailed item list. This will be quite a task to undergo. Maybe I'll start with something (hopefully) simple like z▓*sin(|z|).

The first example of constructing polygons that lie strictly in the black interior region, here for the quadratic set c=-1.125-0.125. Since the black is often disconnected in early refinement steps, usually a set of disjoint polygons is constructed in a similar way as for the outer region. Starting from any 5x5 black square I send rays to all four directions until they encounter gray, step back and mark part of the ray's endpoint as "polygonial".

Since the polygons are constructed to only contain vertical and horizontal lines (I don't like having diagonals when the underlying mathematics in the TSA works on squares parallel to the axis), I could use a simplified version of the even-odd point-in-polygon test algorithm (original paper by Michael Galetzka and Patrick Glauner).

The good thing here is, that the correctness of my implementation can be tested during polygon construction.

For every polygon, the pip-routine checks all image pixels (the polygon lines themselves are conservatively seen as outside). Interior points are marked blue, the polygon is marked yellow. Afterwards the computer checks where in the image the blue lies (so basically a bounding box for blue).

If the blue points only occur within the polygon then the pip-test was implemented correctly for that specific polygon, correct in the sense, that it does not falsely classify an outside pixel as interior. It might not classify all black pixels that lie within it as interior, but that would be type of an error that does not affect the mathematical trustworthiness of the polygon.

Currently some pixels from within are not classified as interior, but that's because I conservatively take the pixel and its 8 neighbours - and only if all nine are within the polygon, then the center pixel is judged as interior.

Attachments:
The first image shows the Julia set with all the identified polygons (yellow) of a decent size.

The second one shows the same set, this time only one polygon is drawn, additionally also all the pixel that the pip-test identified as interior (blue). The intersection rectangle of the red lines contain all the blue pixels - and none is outside the polygon. This test is automated and done purely on the bitmap RGB values for every identified polygon.

The zip file contains the interior polygons as individual point lists, as for the exterior polygon.

Interior and exterior polygons combined could be used as a filter if one is only interested in e.g. the boundary of a Julia set and thus could be used to skip pixels in point sampling.

 

Offline marcm200

  • *
  • Fractal Flamingo
  • ****
  • Posts: 340
« Reply #28 on: June 06, 2019, 09:07:36 AM »
The first trustworthy quadratic Mandelbrot I made. After having the symbolic interval arithmetics parser running, I could fairly simply incorporate claude's suggestion about using intervals for the seed value c as well as for the dynamic z. After 4 days of calculation (with an un-optimized bunch of programs working together) the Mandelbrot in a (low) spatial resolution of 2-9 was finally finished. I'll now discard that routine collection and start again from the beginning to get an optimized, tailor-made program that is hopefully much faster.

Image description:
Axis range -2.+2 in both directions. Directly computed was the region re[-2..+0.8] x im[0..1.3] (white, black, gray), which was then mirrored at the x-axis (blue, dark gray and pink). The region outside [-2..0.8] x [-1.3..1.3] was known to not contain Mandelbrot parts, and was omitted completely from computation (orange, lime), which saved a lot of time.

Offline marcm200

  • *
  • Fractal Flamingo
  • ****
  • Posts: 340
« Reply #29 on: June 24, 2019, 10:24:25 PM »
This quadratic Julia set is quite interesting in how the black interior develops. Usually what I see are some blobs here and there that converge to each other, but here the beginning black shows nicely the underlying spiraling nature.

At a small size of 1024 pixel it looks like a shifted -0.75+0i parabolic set (which interestingly will never develop black squares no matter how deep the refinement goes as the paper mentions), but at 4k the swirling black starts to emerge, which at 16k is nicely developped and shows spirals - a novelty to me.

The image on the upper right shows a zoom into one spiral and the computed exterior (blue) and interior (yellow) polygons which nicely fit into the spirals - something I was very happy to see since during development of the polygon construction routine I only used non-spiraling sets.