4
$\begingroup$

This is my first post here, I have the following basic code for drawing an iterated function system fractal. It works to the third step, but then for the fourth iteration it freezes my old MacBook. I have not included the fourth iteration but it should be clear how to proceed. I reckon it has something to do with the way RegionPlot works. I am looking for an alternative or a way to make it quicker. I have also tried to use DiscretizeRegion to no effect.

ClearAll["Global`*"]
θ = ArcTan[1/2];

sqCreate[{x0_, y0_}, h_] := Rectangle[{x0, y0}, {x0 + h, y0 + h}];
simil[i_, sq_] := Which[
   i == 1, 
   TransformedRegion[sq, AffineTransform[{1/√5 RotationMatrix[-θ],
                                         {0, 1/√5}}]],
   i == 2, 
   TransformedRegion[sq, AffineTransform[{1/√5 RotationMatrix[-θ],
                                         {1/5 + 2/5, 1/√5 - 1/5 + 2/5}}]],
   i == 3, 
   TransformedRegion[sq, AffineTransform[{1/√5 RotationMatrix[-θ],
                                         {2/5 + 2/5, 1/√5 - 1/5 - 1/5}}]],
   i == 4, 
   TransformedRegion[sq, AffineTransform[{1/√5 RotationMatrix[-θ],
                                         {2/5 - 1/5, 1/√5 - 2/5 - 1/5}}]],
   i == 5, 
   TransformedRegion[sq, AffineTransform[{1/√5 RotationMatrix[-θ],
                                         {2/5, 1/√5 - 1/5}}]]
   ];

fullSimil[sq_] := 
  RegionUnion[simil[1, sq], simil[2, sq], simil[3, sq], simil[4, sq], 
   simil[5, sq]];

set0 = sqCreate[{0, 0}, 1];
paint0 = Graphics[{Gray, set0}]

set1 = fullSimil[set0]
paint1 = RegionPlot[set1, PlotRange -> All]

set2 = fullSimil[set1];
paint2 = RegionPlot[set2, PlotRange -> {{0, 1.5}, {-0.5, 1}}, 
  AspectRatio -> Automatic]

set3 = fullSimil[set2];
paint3 = RegionPlot[set3, PlotRange -> {{0, 1.5}, {-1, 1}}, 
  AspectRatio -> Automatic]
$\endgroup$
  • $\begingroup$ Welcome to Mathematica.SE, plebmatician! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour and check the faqs! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! $\endgroup$ – Chris K Jan 4 at 18:52
5
$\begingroup$

This should work a bit faster. The trick is to use BoundaryDiscretizeRegion in each step to convert the internal representation to a single polygonal line. The second ingedient is to use machine precision numbers; these are much faster for numerical computations.

θ = ArcTan[0.5];
A = 1/Sqrt[5] RotationMatrix[-θ];
vecs = Developer`ToPackedArray@N@{{0, 1/Sqrt[5]}, {1/5 + 2/5, 1/Sqrt[5] - 1/5 + 2/5}, {2/5 + 2/5, 1/Sqrt[5] - 1/5 - 1/5}, {2/5 - 1/5, 1/Sqrt[5] - 2/5 - 1/5}, {2/5, 1/Sqrt[5] - 1/5}};

step[sq_] := BoundaryDiscretizeRegion[
  RegionUnion@Table[TransformedRegion[sq, AffineTransform[{A, v}]], {v, vecs}]
  ]

MList = NestList[
  step,
  MeshRegion[N@{{0, 0}, {1, 0}, {1, 1}, {0, 1}}, Polygon[{{1, 2, 3, 4}}]],
  8
  ];

GraphicsGrid[
 Partition[MList, 3],
 ImageSize -> Full
 ]

enter image description here

$\endgroup$
  • $\begingroup$ Thank you very much, amazing how much faster this is. But I don't understand the command "Developer`ToPackedArray", never seen syntax like it before. And what is its purpose, wouldn't the code work the same without it? $\endgroup$ – plebmatician Jan 4 at 19:14
  • 1
    $\begingroup$ Developer`ToPackedArray just converts the array into a so-called packed array. Packed arrays are usually helpful for numerical performance. But here Developer`ToPackedArray makes probably no difference. $\endgroup$ – Henrik Schumacher Jan 4 at 19:29
2
$\begingroup$

With a little refactoring of your code, and with the use of some undocumented functionality:

θ = ArcTan[1/2];
sqCreate[{x0_, y0_}, h_] := Polygon[{{x0, y0}, {x0 + h, y0}, {x0 + h, y0 + h}, {x0, y0 + h}}]

fullSimil[sq_Polygon] := Map[First, Graphics`PolygonUtils`PolygonCombine[
                      Table[Polygon[AffineTransform[{RotationMatrix[-θ]/Sqrt[5], pos}] @@ sq],
                            {pos, {{0, 1/Sqrt[5]}, {3/5, 1/Sqrt[5] + 1/5}, {4/5, 1/Sqrt[5] - 2/5},
                                   {1/5, 1/Sqrt[5] - 3/5}, {2/5, 1/Sqrt[5] - 1/5}}}]]]

Partition[BoundaryDiscretizeGraphics /@ NestList[fullSimil, sqCreate[{0, 0}, 1], 8],
          3] // GraphicsGrid

fractal meshes

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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