Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

OP edit: This is a Mathematica-specific question about an approach it attempted for a fractal visualization problem described HERE. I'm using the Mandelbrot set there and here as an example, but the solutions apply to a broader fractal problem.

I need to apply an iterative function z to a list.

z1[n_, c_] := If[n > 0, z1[n - 1, c]^2 + c, c];
c = {-1, -.5, 0, .5, 1};
z1[7, c]

(* {0, -0.350234, 0, 12005.5, 44127887745906175987802} *)

It happens that if Abs[z[n, c]] > 2 at n, then all subsequent evaluations are also greater than 2. I only need to know the values of c for which the z series satisfies Abs[z] < 2.

I tried to make my calculation more efficient by defining my function like this:

z2[n_, c_] := NestWhile[(#^2 + c) &, c, Abs[#] <= 2 &, 1, n];

This works great on single values of c:

z2[12, 1]

(* 5 *)

but it just returns the List back to me if c is a List:

z2[12, c]

(* {-1, -0.5, 0, 0.5, 1} *)

How can I fix this? I don't want to use a loop because I need to process long lists to high iterations values. In the example below, a lot of time is wasted iterating values that already satisfy Abs[z] < 2.

stepSize = .0001;
iter = 46;
c = Flatten[Table[x + I*y, {x, -1.42, -1.39, stepSize}, {y, -.005, .025, stepSize}]];
p = z1[iter, c];
r = Pick[c, Abs[#] <= 2 & /@ p];
ListPlot[Transpose[{Re[r], Im[r]}]]

Fractal

share|improve this question
1  
Related: (103673) (yes, your own question; I'm linking it for others.) – Mr.Wizard Jan 16 at 22:43
    
Ah, @Mr.Wizard beat me to linking that. The code I tried for z2[] worked perfectly for a ContourPlot but doesn't work on a List. – Jerry Guern Jan 16 at 22:55
    
The If requires a scalar argument. z2[12,#]& /@ c should work. – mikado Jan 16 at 23:37
    
This sounds like job for CUDA, about which I am largely ignorant. Then plotting via ListContourPlot[]. – Eric Towers Jan 17 at 5:24

You can do it probably most efficiently in compiled code, if you're not too concerned about precision. Here you can use the listability of compiled functions over tensor arguments.

Your function is basically the Mandelbrot iteration:

mandelbrot = Compile[{{c, _Complex, 0}, {d, _Integer, 0}},
   Block[{i = 0, z = c}, While[Abs[z] < 2.0 && i < d, z = z*z + c; i++]; i],
   RuntimeOptions -> "Speed",
   RuntimeAttributes -> Listable, Parallelization -> True
  ];

AbsoluteTiming[
 result = mandelbrot[
    Outer[Complex, Range[-1.42, -1.39, .0001], Range[-.005, .025, .0001]],
    80
  ];
 ] (* -> 0.219 seconds *)

ArrayPlot[result, Frame -> False, PlotRangePadding -> None, PixelConstrained -> {1, 1}]

plot of Mandelbrot set

If precision is very important, you will not be able to use compiled code, so the best you can do is probably just have to give your function the Listable attribute. This is not particularly efficient; it just uses Thread when the function is called with lists in the arguments.

share|improve this answer
    
Your comments intrigue me! Why can't I use compiled code if precision is important? Precision is very important in what I'm doing, because I need to see the filament structure of my fractals, not the color show around them. That's why my plots are monochrome. – Jerry Guern Jan 17 at 5:40
    
When you use Compile, then the numbers will always be machine precision numbers like double in C/C++. Mathematica is able to calculate with arbitrary precision numbers, but as soon as you compile something, you are stuck with the usual machine type numbers. – halirutan Jan 17 at 11:24
    
@JerryGuern the "colors" record, for each point, the number of iterations performed (i.e., the number required for the value to reach the basin of attraction). The result is a contour plot of exactly that. I thought that's what you wanted? halirutan's comment explains why compiled code cannot be used for high precision calculations and I don't have anything to add to that. – Oleksandr R. Jan 17 at 15:41
    
@JerryGuern as we discussed on your last question, you will have a limited number of digits of precision left after a certain number of iterations. Starting from machine precision, you will have 11 digits left after 100 iterations; 6 after 200; and 2 after 400. So it depends on what your precision requirements are and how many iterations you want to do (and somewhat on the starting point) as to whether machine precision will be enough or not. – Oleksandr R. Jan 17 at 16:00
2  
@OleksandrR. About the "colors". I don't ultimately want to display the contours that lay outside the mandelbrot set, I was to display just the set itself but with the filaments exaggerated in thickness just enough to keep them visible and the displayed set visibly connected. That's what led to my playing with ContourPlot. If I could display the z1[n,c]==2 contour and increase n so that it's indistinguishable from the n->Inf contour at current resolution, that would give me what I want. Not a great idea, I know, but it motivated my questions about ContourPlot. – Jerry Guern Jan 17 at 20:43

Again I will propose solving the problem by a recursive function.

helper[c0_, c_] :=
  Module[{r = (Abs[#] < 2 & /@ c)}, Pick[#, r] & /@ {c0, c^2 + c0}]
z[n_Integer?Positive, c0_List] := z[n - 1, c0, c0^2 + c0]
z[0, c0_, _] := c0
z[n_, c0_, c_] := z[n - 1, Sequence @@ helper[c0, c]]

With[{stepSize = .0001, n = 46},
  Module[{c0, p},
    c0 = 
      Flatten[
        Table[x + I*y, 
          {x, -1.42, -1.39, stepSize}, {y, -.005, .025, stepSize}]]; 
    p = z[n, c0];
    ListPlot[Transpose[{Re[p], Im[p]}]]]]

plot

Note: because z is tail recursive, Mathematica eliminates the recursion and evaluates it iteratively. So this should be considered an iterative solution as requested.

share|improve this answer
    
Ah. I always learn a lot from your answers, @m_goldberg. But why did you put mandelbrot in the title? I'm using the mandelbrot iterator as an example because it's simpler than my iterator, but the techniques being discussed aren't limited to the mandelbrot set. – Jerry Guern Jan 17 at 5:42
1  
@JerryGuern. I put "Mandelbrot" in the title because I think this question and both answers will be of interest to people looking for Mathematica generators of the Mandelbrot set. I did it to make it easy for them to find the question in searches. I saw nothing in your question that indicated to me you were really interested in iteration in the abstract. – m_goldberg Jan 17 at 9:09
1  
@JerryGuern I believe it may be useful to you to search on this site using "mandelbrot" as a search term. The subject of fast iterations has been discussed previously. Perhaps I should mention that, for me at least, with your multiple very similar questions about the subject it is becoming a little difficult to discern what you really want. – Oleksandr R. Jan 17 at 17:18
    
@OleksandrR. I've done that. I am working on a fractal problem with similarities to the mandelbrot set. But each of my similar recent questions posted here has been about a distinct specific MMa issue I encountered while working on the fractal problem. I've posted math-specific questions on the math SE site. I'm trying hard to stick to the rules and keep my questions distinct, specific, and on the appropriate site. I've had to do this because SE is a Q-A site not a discussion site. – Jerry Guern Jan 17 at 20:07
2  
@JerryGuern your recent questions are certainly very concrete, but it seems to me that answers risk missing the point without some hint to your overall purpose. Your explanation here makes sense and I understand now why you need high precision, so thank you for that. My suggestion would be not to be extremely concerned with "the rules" but rather simply to ensure that your questions cannot too easily be misinterpreted in such a way as to conflict with your real goals. Being as specific as possible is an important part of that, IMO. – Oleksandr R. Jan 17 at 21:08

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.