3
$\begingroup$

Consider the following code:

c0 := -1.24; 
b := 0.125; 
hh[{x_, y_}] := {x^2 + c0 - b*y, x}; 
hh2[{x_, y_}] := hh[hh[{x, y}]]; 
dhh := Transpose[{D[hh[{x, y}], x], D[hh[{x, y}], y]}]; 
dhh2 := Transpose[{D[hh[hh[{x, y}]], x], D[hh[hh[{x, y}]], y]}]; 
fixpoints = NSolve[hh[{x, y}] == {x, y}, {x, y}]; 
period2cycle = NSolve[hh[hh[{x, y}]] == {x, y}, {x, y}]; 
center = {x, y} /. fixpoints[[1]]; 
centerOther = {x, y} /. fixpoints[[2]]; 
pushvector = Eigensystem[dhh /. fixpoints[[1]]][[2,1]]; 
pushvectorOther = Eigensystem[dhh /. fixpoints[[2]]][[2,1]]; 
sink = {x, y} /. period2cycle[[1]]; 
pushvector2 = Eigensystem[dhh /. fixpoints[[2]]][[2,1]]; 
maxIter = 63; 
maxXXSize = 600; 
maxYYSize = 600; 
smalldelta = 0.05/maxXXSize; 
newfun[{{x_, y_}, n_}] := {hh[{x, y}], n + 1}; 
newtest[{{x_, y_}, n_}] := Abs[x] < 10000 && 
    Norm[{x, y} - sink] > 0.001; 
newcolorassign[{{x_, y_}, n_}] := Which[Abs[x] < 3, {0, 0, 0}, 
    Mod[Floor[Arg[x]/2^16], 2] > 0, {0.5, 0.5, 0.5}, True, 
    {1, 1, 1}]; 

With this setup one can execute the following:

Graphics[{Raster[Table[newcolorassign[NestWhile[newfun, 
       {center - pushvector*smalldelta*I*(k - maxXXSize/2 + 
           I*(j - maxYYSize/2)), 0}, newtest, 1, maxIter]], 
     {k, 1, maxXXSize}, {j, 1, maxYYSize}]], 
   PointSize[12/maxXXSize], White, 
   Point[{(maxXXSize - 1)/2, (maxYYSize - 1)/2}], 
   PointSize[4/maxXXSize], Black, 
   Point[{(maxXXSize - 1)/2, (maxYYSize - 1)/2}]}]

In a notebook this will produce a picture but because you specify the parameters at the beginning, to avoid confusion you should start a new notebook for each parameter pair (b,c0). Now to get around this problem I tried to generalise the code to the following.

hh1[{x_, y_}, b_, c_] := {x^2 + c - b*y, x}; 
hhhh[{x_, y_}, b_, c_] := hh1[hh1[{x, y}, b, c], b, c]; 
dhhh[b_, c_] := Transpose[{D[hh1[{x, y}, b, c], x], 
     D[hh1[{x, y}, b, c], y]}]; 
dhhh2[b_, c_] := Transpose[{D[hhhh[{x, y}, b, c], x], 
     D[hhhh[{x, y}, b, c], y]}]; 
fixpts[b_, c_] := NSolve[hh1[{x, y}, b, c] == {x, y}, {x, y}]
per2cycle[b_, c_] := NSolve[hhhh[{x, y}, b, c] == {x, y}, {x, y}]
centre[b_, c_] := {x, y} /. fixpts[b, c][[1]]; 
centreOther[b_, c_] := {x, y} /. fixpts[b, c][[2]]; 
pushvect[b_, c_] := Eigensystem[dhhh[b, c] /. 
      fixpts[b, c][[1]]][[2,1]]; 
pushvectOther[b_, c_] := 
   Eigensystem[dhhh[b, c] /. fixpts[b, c][[2]]][[2,1]]; 
snk[b_, c_] := {x, y} /. per2cycle[b, c][[1]]; 
pushvect2[b_, c_] := Eigensystem[dhhh[b, c] /. 
      fixpts[b, c][[2]]][[2,1]]; 
maxIter = 63; 
maxXXSize = 600; 
maxYYSize = 600; 
smalldelta = 0.05/maxXXSize; 
newfun[{{x_, y_}, n_}, b_, c_] := {hh1[{x, y}, b, c], n + 1}; 
newtest[{{x_, y_}, n_}, b_, c_] := Abs[x] < 10000 && 
    Norm[{x, y} - snk[b, c]] > 0.001; 
newcolorassign[{{x_, y_}, n_}] := Which[Abs[x] < 3, {0, 0, 0}, 
    Mod[Floor[Arg[x]/2^16], 2] > 0, {0.5, 0.5, 0.5}, True, 
    {1, 1, 1}]; 
unstablepic[b_, c_] := Graphics[
    {Raster[Table[newcolorassign[NestWhile[newfun[#1, b, c] & , 
         {centre[b, c] - pushvect[b, c]*smalldelta*I*
            (k - maxXXSize/2 + I*(j - maxYYSize/2)), 0}, 
         newtest[#1, b, c] & , 1, maxIter]], {k, 1, maxXXSize}, 
       {j, 1, maxYYSize}]], PointSize[12/maxXXSize], White, 
     Point[{(maxXXSize - 1)/2, (maxYYSize - 1)/2}], 
     PointSize[4/maxXXSize], Black, 
     Point[{(maxXXSize - 1)/2, (maxYYSize - 1)/2}]}]; 

Now unstablepic[0.125,-1.24] produces the same output as the original code but is significantly slower. To see this more quickly, one can take maxXXSize = 60 and maxYYSize = 60 instead of 600. I'm not sure why the generalised code is taking so much longer than the original. Is there a way to make this quicker?

$\endgroup$
6
$\begingroup$

It looks to me like you're recalculating a lot of stuff that should be constant (unless you're doing dangerous things with global variables).

Here's a way to not do that, mostly by changing newTest and

newtest[{{x_, y_}, n_}, b_, c_, snk_] :=
  Abs[x] < 10000 && Norm[{x, y} - snk] > 0.001;
unstablePts[b_, c_] :=
 With[
  {
   centre = centre[b, c], 
   pushvect = pushvect[b, c],
   sink = snk[b, c]
   },
  Table[
   NestWhile[
    newfun[#1, b, c] &,
    {
     centre - pushvect*
       smalldelta*I*(k - maxXXSize/2 + I*(j - maxYYSize/2)), 0},
    newtest[#1, b, c, sink] &,
    1,
    maxIter
    ],
   {k, 1, maxXXSize},
   {j, 1, maxYYSize}
   ]
  ]

unstablepic[b_, c_] :=
  Graphics[{
    Raster[Map[newcolorassign, unstablePts[b, c], {2}]],
    PointSize[12/maxXXSize], White, 
    Point[{(maxXXSize - 1)/2, (maxYYSize - 1)/2}], PointSize[4/maxXXSize], 
    Black, Point[{(maxXXSize - 1)/2, (maxYYSize - 1)/2}]
    }];

you can be a bit more efficient, though, and Compile part of this, (don't worry too much about the "Formal" variables, I'm just using them in case you accidentally defined x, y, b, or c somewhere)

hh1C =
  With[{code = hh1[{\[FormalX], \[FormalY]}, \[FormalB], \[FormalC]]},
   Compile[{{\[FormalX], _Complex}, {\[FormalY], _Complex}, {\[FormalB], \
_Real}, {\[FormalC], _Real}}, 
    code,
    CompilationTarget -> "C"
    ]
   ];
newfunC[{{x_, y_}, n_}, b_, c_] := {hh1C[x, y, b, c], n + 1};

the speed up is relatively minor, but will increase with image size

Block[{
   maxIter = 63,
   maxXXSize = 500,
   maxYYSize = 500
   },
  unstablepic[0.125, -1.24];
  ] // AbsoluteTiming

{20.7144, Null}

Block[{
    newfun = newfunC,
    maxIter = 63,
    maxXXSize = 500,
    maxYYSize = 500
    },
   pic = unstablepic[0.125, -1.24];
   ] // AbsoluteTiming // First

18.9731

pic

enter image description here

$\endgroup$
0

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.