{"id":773,"date":"2013-09-30T12:00:01","date_gmt":"2013-09-30T17:00:01","guid":{"rendered":"https:\/\/blogs.mathworks.com\/cleve\/?p=773"},"modified":"2018-05-25T14:52:52","modified_gmt":"2018-05-25T19:52:52","slug":"iterated-powers-fractal","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/cleve\/2013\/09\/30\/iterated-powers-fractal\/","title":{"rendered":"Iterated Powers Fractal"},"content":{"rendered":"<div class=\"content\">\r\n\r\n<!--introduction-->Iterating powers of complex numbers can produce cycles. Colormaps of the lengths of these cycles produce fractal images.\r\n\r\n<img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/cleve\/view1.png\" alt=\"\" hspace=\"5\" vspace=\"5\" \/>\r\n\r\n<!--\/introduction-->\r\n<h3>Contents<\/h3>\r\n<div>\r\n<ul>\r\n \t<li><a href=\"#9e2cd352-1fcb-472a-a39e-59cb8f5684cd\">UWO Poster<\/a><\/li>\r\n \t<li><a href=\"#c531e9c3-d38b-4c7f-8efa-57d33e00d776\">The Iteration<\/a><\/li>\r\n \t<li><a href=\"#ace40bb4-3d04-4eac-aa89-775718e95400\">An Example<\/a><\/li>\r\n \t<li><a href=\"#6c66a954-1c19-4cf5-ac13-170e196d27dd\">Cycle Length<\/a><\/li>\r\n \t<li><a href=\"#2f90b507-a178-46b0-900d-37fbab16b7f6\">Colormaps<\/a><\/li>\r\n \t<li><a href=\"#ee620a04-372b-480e-ab19-4a163221b044\">Generate an Image<\/a><\/li>\r\n \t<li><a href=\"#09864b3c-57a2-4ae3-a4aa-e4cd056d3c76\">First Image<\/a><\/li>\r\n \t<li><a href=\"#955bdd49-6504-4313-b6d1-4ec44c60ebdd\">Second Image<\/a><\/li>\r\n \t<li><a href=\"#cf2b7945-b513-4cf0-b5b1-b254b0fe0705\">Third Image<\/a><\/li>\r\n \t<li><a href=\"#177f195a-29f2-4af6-b4ff-aa89855cb3d4\">Reference<\/a><\/li>\r\n<\/ul>\r\n<\/div>\r\n<h4>UWO Poster<a name=\"9e2cd352-1fcb-472a-a39e-59cb8f5684cd\"><\/a><\/h4>\r\nIn <a href=\"https:\/\/blogs.mathworks.com\/cleve\/2013\/09\/02\/the-lambert-w-function\">my blog<\/a> a month ago I mentioned <a title=\"http:\/\/www.orcca.on.ca\/LambertW\/ (link no longer working)\">a poster<\/a> produced at the University of Western Ontario about the Lambert W function. And in <a href=\"https:\/\/blogs.mathworks.com\/cleve\/2013\/09\/16\/tower-of-powers\/\">my blog<\/a> two weeks ago about the tower of powers\r\n\r\n$$ z^{z^{z^z}} $$\r\n\r\nI said that if the base $z$ is less that $\\mbox{e}^{-\\mbox{e}}$, which is about 0.066, then the iteration evaluating the tower settles into a stable two-cycle. Both of these posts lead to today's blog where $z$ moves in the complex plane.\r\n<h4>The Iteration<a name=\"c531e9c3-d38b-4c7f-8efa-57d33e00d776\"><\/a><\/h4>\r\nThe base of the iteration is a real or complex number $z$ that stays fixed throughout the iteration. Start with $y_0 = 1$ and simply iterate\r\n\r\n$$ y_{k+1} = z^{y_k} $$\r\n\r\nIf $z$ is too large, then $y_k$ will quickly blow up. As we saw in my last blog, this happens for real $z$ larger than $\\mbox{e}^{1\/\\mbox{e}}$, which is about $1.444$.\r\n\r\nFor some $z$, the iterates $y_k$ approach a fixed point. An example is $z = \\sqrt{2}$. The $y_k$ converge to $2$. We regard fixed points as cycles of length 1.\r\n<h4>An Example<a name=\"ace40bb4-3d04-4eac-aa89-775718e95400\"><\/a><\/h4>\r\nHere is an example of a complex $z$ that produces a cycle of length 7, $z = 2.5+4i$.\r\n<pre class=\"codeinput\">   format <span class=\"string\">short<\/span>\r\n   z = 2.5+4i\r\n   y = 1\r\n   <span class=\"keyword\">for<\/span> k = 1:7\r\n      y = z^y;\r\n      disp(y)\r\n   <span class=\"keyword\">end<\/span>\r\n<\/pre>\r\n<pre class=\"codeoutput\">z =\r\n   2.5000 + 4.0000i\r\ny =\r\n     1\r\n   2.5000 + 4.0000i\r\n  -0.6503 + 0.5363i\r\n   0.2087 + 0.0366i\r\n   1.2845 + 0.3528i\r\n  -1.4011 + 4.9361i\r\n   7.6883e-04 - 3.4362e-05i\r\n   1.0012 + 0.0007i\r\n<\/pre>\r\nAfter 7 steps, we're not quite back at <tt>y = 1<\/tt>, but we're pretty close. Repeating the iteration a few dozen times brings us to a stable cycle of length 7 near this initial transient. Here is a picture.\r\n\r\n<img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/cleve\/cycleplot.png\" alt=\"\" hspace=\"5\" vspace=\"5\" \/>\r\n<h4>Cycle Length<a name=\"6c66a954-1c19-4cf5-ac13-170e196d27dd\"><\/a><\/h4>\r\nFinding cycle length involves saving the iterates on a stack and comparing each new iterate with the previous ones. The internal parameters of the function are the tolerance in the comparison, the stack size, and the maximum number of iterations. The cycle length is declared to be infinite if either the iterates go off to infinity or no match is found before the iteration limit is reached.\r\n\r\nHere is the function.\r\n<pre class=\"codeinput\">   type <span class=\"string\">cycle<\/span>\r\n<\/pre>\r\n<pre class=\"codeoutput\">function c = cycle(z)\r\n% CYCLE  Cycle length of iterated power\r\n% c = cycle(z)\r\n   tol = 1.e-4;   % Tolerance\r\n   N = 40;        % Height of stack\r\n   M = 2000;      % Maximum number of iterations\r\n   S = zeros(N,1);   % Stack\r\n   S(:) = NaN;\r\n   y = 1;\r\n   k = 0;\r\n   c = Inf;\r\n   while k &lt; M\r\n      y = z^y;\r\n      k = k+1;\r\n      if isinf(y)\r\n         break\r\n      elseif any(abs(y-S) &lt; tol)\r\n         c = min(find(abs(y-S) &lt; tol));\r\n         break\r\n      end\r\n      S(2:N) = S(1:N-1);\r\n      S(1) = y;\r\n   end\r\nend\r\n\r\n<\/pre>\r\nLet's test the function on a few cases where we already know the answer.\r\n<pre class=\"codeinput\">   cycle(sqrt(2))  <span class=\"comment\">% Should be 1<\/span>\r\n   cycle(1.5)      <span class=\"comment\">% Should be Inf<\/span>\r\n   cycle(.05)      <span class=\"comment\">% Should be 2<\/span>\r\n   cycle(2.5+4i)   <span class=\"comment\">% Should be 7<\/span>\r\n<\/pre>\r\n<pre class=\"codeoutput\">ans =\r\n     1\r\nans =\r\n   Inf\r\nans =\r\n     2\r\nans =\r\n     7\r\n<\/pre>\r\n<h4>Colormaps<a name=\"2f90b507-a178-46b0-900d-37fbab16b7f6\"><\/a><\/h4>\r\nGood colormaps are an important aspect of generating good fractal images. I am using the colormaps from the <a href=\"https:\/\/www.mathworks.com\/content\/dam\/mathworks\/mathworks-dot-com\/moler\/exm\/chapters\/mandelbrot.pdf\">Mandelbrot chapter<\/a> of <a href=\"https:\/\/www.mathworks.com\/moler\/exm\/\">Experiments with MATLAB<\/a>. These take small chunks of the colormaps from MATLAB itself and repeat them several times. The maps are called <tt>jets<\/tt>, <tt>hots<\/tt>, and <tt>cmyk<\/tt> and are included in the <a href=\"https:\/\/www.mathworks.com\/moler\/exm\/chapters.html\">mandelbrot.m code<\/a>. Here is one of them.\r\n<pre class=\"codeinput\">   type <span class=\"string\">jets<\/span>\r\n<\/pre>\r\n<pre class=\"codeoutput\">function cmap = jets(m)\r\n   c = jet(16);\r\n   e = ones(m\/16,1);\r\n   cmap = kron(e,c);\r\nend\r\n\r\n<\/pre>\r\n<h4>Generate an Image<a name=\"ee620a04-372b-480e-ab19-4a163221b044\"><\/a><\/h4>\r\nTo generate a fractal image, put a grid over a region in the complex region, compute the cycle length for each grid point, and then use the MATLAB <tt>image<\/tt> function and one of the colormaps to produce a pseudocolor image. (This is an <i>embarrassingly parallel<\/i> task and has for many years been used to demonstrate parallel computers. Load balancing is a challenge. But I digress ...)\r\n<h4>First Image<a name=\"09864b3c-57a2-4ae3-a4aa-e4cd056d3c76\"><\/a><\/h4>\r\nThe image at the top of this post shows a large chunk of the complex plane. With the <tt>cmyk<\/tt> colormap, the infinite cycle lengths come out black. I can see that with a <tt>spy<\/tt> plot on my machine where I have saved the data.\r\n<pre class=\"codeinput\">   load <span class=\"string\">view1<\/span>\r\n   spy(isinf(u))\r\n   set(gca,<span class=\"string\">'ydir'<\/span>,<span class=\"string\">'n'<\/span>)\r\n<\/pre>\r\n<img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/cleve\/ipf_01.png\" alt=\"\" hspace=\"5\" vspace=\"5\" \/>\r\n\r\nWhile I am it, let's check the density of <tt>Inf's.<\/tt>\r\n<pre class=\"codeinput\">   nnz(isinf(u))\/prod(size(u))\r\n<\/pre>\r\n<pre class=\"codeoutput\">ans =\r\n    0.2602\r\n<\/pre>\r\nOver a quarter of the grid points have infinite cycle lengths.\r\n<h4>Second Image<a name=\"955bdd49-6504-4313-b6d1-4ec44c60ebdd\"><\/a><\/h4>\r\nZoom in on a smaller region in the plane, and use the <tt>jets<\/tt> colormap with its columns flipped. This is typical of the structure we see everywhere in this fractal, at all magnification levels. They are certainly not as rich as the Mandelbrot set, but I still them interesting.\r\n\r\n<img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/cleve\/view2.png\" alt=\"\" hspace=\"5\" vspace=\"5\" \/>\r\n<h4>Third Image<a name=\"cf2b7945-b513-4cf0-b5b1-b254b0fe0705\"><\/a><\/h4>\r\nZoom in farther, and use the <tt>hots<\/tt> colormap. The large brown region in the northwest corner is not infinite cycle length. It happens to be <tt>cycle == 16<\/tt>.\r\n\r\n<img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/cleve\/view3.png\" alt=\"\" hspace=\"5\" vspace=\"5\" \/>\r\n\r\nHere is the code that generates this image. It takes about 5 minutes to run on my laptop, which is a couple of years old.\r\n<pre class=\"codeinput\">   type <span class=\"string\">view3<\/span>\r\n<\/pre>\r\n<pre class=\"codeoutput\">x0 = 2.3825;\r\ny0 = 3.7227;\r\nd = 1.e-5;\r\nnhalf = 200;\r\nv = d*(-nhalf:nhalf);\r\n[x,y] = meshgrid(x0+v,y0+v);\r\nz = x+i*y;\r\nu = zeros(size(z));\r\n\r\nn = size(z,1)\r\ntic\r\nfor k = 1:n\r\n   for j = 1:n\r\n      u(k,j) = cycle(z(k,j));\r\n   end\r\n   k\r\nend\r\ntoc\r\n\r\nimage(u)\r\naxis square\r\ncolormap(hots(128))\r\nset(gca,'ydir','n')\r\nset(gca,'xtick',[1,n\/2,n])\r\nset(gca,'ytick',[1,n\/2,n])\r\nset(gca,'xticklabel',{min(x(:)),num2str(x0),max(x(:))})\r\nset(gca,'yticklabel',{min(y(:)),num2str(y0),max(y(:))})\r\nxlabel('x')\r\nylabel('y')\r\n\r\n<\/pre>\r\n<h4>Reference<a name=\"177f195a-29f2-4af6-b4ff-aa89855cb3d4\"><\/a><\/h4>\r\nThe Lambert W Function, Poster, Ontario Research Centre for Computer Algebra, http:\/\/www.orcca.on.ca\/LambertW\/\r\n\r\n<script>\/\/ <![CDATA[\r\nfunction grabCode_31a7493b75a347eabf4323b8a302609f() {\r\n        \/\/ Remember the title so we can use it in the new page\r\n        title = document.title;\r\n\r\n        \/\/ Break up these strings so that their presence\r\n        \/\/ in the Javascript doesn't mess up the search for\r\n        \/\/ the MATLAB code.\r\n        t1='31a7493b75a347eabf4323b8a302609f ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' 31a7493b75a347eabf4323b8a302609f';\r\n    \r\n        b=document.getElementsByTagName('body')[0];\r\n        i1=b.innerHTML.indexOf(t1)+t1.length;\r\n        i2=b.innerHTML.indexOf(t2);\r\n \r\n        code_string = b.innerHTML.substring(i1, i2);\r\n        code_string = code_string.replace(\/REPLACE_WITH_DASH_DASH\/g,'--');\r\n\r\n        \/\/ Use \/x3C\/g instead of the less-than character to avoid errors \r\n        \/\/ in the XML parser.\r\n        \/\/ Use '\\x26#60;' instead of '<' so that the XML parser\r\n        \/\/ doesn't go ahead and substitute the less-than character. \r\n        code_string = code_string.replace(\/\\x3C\/g, '\\x26#60;');\r\n\r\n        copyright = 'Copyright 2013 The MathWorks, Inc.';\r\n\r\n        w = window.open();\r\n        d = w.document;\r\n        d.write('\r\n\r\n\r\n\r\n<pre>\\n');\r\n        d.write(code_string);\r\n\r\n        \/\/ Add copyright line at the bottom if specified.\r\n        if (copyright.length > 0) {\r\n            d.writeln('');\r\n            d.writeln('%%');\r\n            if (copyright.length > 0) {\r\n                d.writeln('% _' + copyright + '_');\r\n            }\r\n        }\r\n\r\n        d.write('<\/pre>\r\n\r\n\r\n\r\n\r\n\\n');\r\n\r\n        d.title = title + ' (MATLAB code)';\r\n        d.close();\r\n    }\r\n\/\/ ]]><\/script>\r\n<p style=\"text-align: right; font-size: xx-small; font-weight: lighter; font-style: italic; color: gray;\"><a><span style=\"font-size: x-small; font-style: italic;\">Get\r\nthe MATLAB code<noscript>(requires JavaScript)<\/noscript><\/span><\/a><\/p>\r\nPublished with MATLAB\u00ae R2013b\r\n<p class=\"footer\">Published with MATLAB\u00ae R2013b<\/p>\r\n\r\n<\/div>\r\n<!--\r\n31a7493b75a347eabf4323b8a302609f ##### SOURCE BEGIN #####\r\n%% Iterated Powers Fractal\r\n% Iterating powers of complex numbers can produce cycles.\r\n% Colormaps of the lengths of these cycles produce fractal images.\r\n%\r\n% <<view1.png>>\r\n%\r\n\r\n%% UWO Poster\r\n% In <https:\/\/blogs.mathworks.com\/cleve\/2013\/09\/02\/the-lambert-w-function % my blog> a month ago I mentioned <http:\/\/www.orcca.on.ca\/LambertW\/ a poster>\r\n% produced at the University of Western Ontario about the Lambert W function.\r\n% And in\r\n% <https:\/\/blogs.mathworks.com\/cleve\/2013\/09\/16\/tower-of-powers\/ % my blog> two weeks ago about the tower of powers\r\n%\r\n% $$ z^{z^{z^z}} $$\r\n%\r\n% I said that if the base $z$ is less that $\\mbox{e}^{-\\mbox{e}}$,\r\n% which is about 0.066, then the iteration evaluating the tower settles\r\n% into a stable two-cycle.  Both of these posts lead to today's blog\r\n% where $z$ moves in the complex plane.\r\n\r\n%% The Iteration\r\n% The base of the iteration is a real or complex number $z$ that stays\r\n% fixed throughout the iteration.  Start with $y_0 = 1$ and simply iterate\r\n%\r\n% $$ y_{k+1} = z^{y_k} $$\r\n%\r\n% If $z$ is too large, then $y_k$ will quickly blow up.  As we saw in my\r\n% last blog, this happens for real $z$ larger than $\\mbox{e}^{1\/\\mbox{e}}$,\r\n% which is about $1.444$.\r\n\r\n%%\r\n% For some $z$, the iterates $y_k$ approach a fixed point.\r\n% An example is $z = \\sqrt{2}$.  The $y_k$ converge to $2$.\r\n% We regard fixed points as cycles of length 1.\r\n\r\n%% An Example\r\n% Here is an example of a complex $z$ that produces a cycle of length 7,\r\n% $z = 2.5+4i$.\r\n%\r\nformat short\r\nz = 2.5+4i\r\ny = 1\r\nfor k = 1:7\r\ny = z^y;\r\ndisp(y)\r\nend\r\n\r\n%%\r\n% After 7 steps, we're not quite back at |y = 1|, but we're pretty close.\r\n% Repeating the iteration a few dozen times brings us to a stable cycle\r\n% of length 7 near this initial transient.  Here is a picture.\r\n%\r\n% <<cycleplot.png>>\r\n%\r\n\r\n%% Cycle Length\r\n% Finding cycle length involves saving the iterates on a stack and\r\n% comparing each new iterate with the previous ones.  The internal parameters\r\n% of the function are the tolerance in the comparison, the stack\r\n% size, and the maximum number of iterations.  The cycle length is\r\n% declared to be infinite if either the iterates go off to infinity or\r\n% no match is found before the iteration limit is reached.\r\n\r\n%%\r\n% Here is the function.\r\n%\r\ntype cycle\r\n\r\n%%\r\n% Let's test the function on a few cases where we already know the answer.\r\n%\r\ncycle(sqrt(2))  % Should be 1\r\ncycle(1.5)      % Should be Inf\r\ncycle(.05)      % Should be 2\r\ncycle(2.5+4i)   % Should be 7\r\n\r\n%% Colormaps\r\n% Good colormaps are an important aspect of generating good fractal images.\r\n% I am using the colormaps from the\r\n% <https:\/\/www.mathworks.com\/moler\/exm\/chapters\/mandelbrot.pdf % Mandelbrot chapter> of <https:\/\/www.mathworks.com\/moler\/exm\/ % Experiments with MATLAB>.  These take small chunks of the colormaps from\r\n% MATLAB itself and repeat them several times.  The maps are called\r\n% |jets|, |hots|, and |cmyk| and are included in the\r\n% <https:\/\/www.mathworks.com\/moler\/exm\/exmfilelist.html mandelbrot.m code>.\r\n% Here is one of them.\r\n%\r\ntype jets\r\n\r\n%% Generate an Image\r\n% To generate a fractal image, put a grid over a region in the complex\r\n% region, compute the cycle length for each grid point, and then use\r\n% the MATLAB |image| function and one of the colormaps to produce a\r\n% pseudocolor image.  (This is an _embarrassingly parallel_ task and has\r\n% for many years been used to demonstrate parallel computers.  Load\r\n% balancing is a challenge.  But I digress ...)\r\n\r\n%% First Image\r\n% The image at the top of this post shows a large chunk of the complex plane.\r\n% With the |cmyk| colormap, the infinite cycle lengths come out black.\r\n% I can see that with a |spy| plot on my machine where I have saved the data.\r\n%\r\nload view1\r\nspy(isinf(u))\r\nset(gca,'ydir','n')\r\n\r\n%%\r\n% While I am it, let's check the density of |Inf's.|\r\n%\r\nnnz(isinf(u))\/prod(size(u))\r\n\r\n%%\r\n% Over a quarter of the grid points have infinite cycle lengths.\r\n%\r\n\r\n%% Second Image\r\n% Zoom in on a smaller region in the plane, and use the |jets| colormap\r\n% with its columns flipped.\r\n% This is typical of the structure we see everywhere in this fractal,\r\n% at all magnification levels.\r\n% They are certainly not as rich as the Mandelbrot set, but I still them\r\n% interesting.\r\n%\r\n% <<view2.png>>\r\n%\r\n\r\n%% Third Image\r\n% Zoom in farther, and use the |hots| colormap.  The large brown region\r\n% in the northwest corner is not infinite cycle length.  It happens to be\r\n% |cycle == 16|.\r\n%\r\n% <<view3.png>>\r\n%\r\n\r\n%%\r\n% Here is the code that generates this image.\r\n% It takes about 5 minutes to run on my laptop, which is a couple\r\n% of years old.\r\n%\r\ntype view3\r\n\r\n%% Reference\r\n% The Lambert W Function, Poster, Ontario Research Centre for Computer Algebra,\r\n% <http:\/\/www.orcca.on.ca\/LambertW\/ http:\/\/www.orcca.on.ca\/LambertW>\r\n\r\n##### SOURCE END ##### 31a7493b75a347eabf4323b8a302609f\r\n-->","protected":false},"excerpt":{"rendered":"<div class=\"overview-image\"><img decoding=\"async\"  class=\"img-responsive\" src=\"https:\/\/blogs.mathworks.com\/images\/cleve\/view1.png\" onError=\"this.style.display ='none';\" \/><\/div>... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/cleve\/2013\/09\/30\/iterated-powers-fractal\/\">read more >><\/a><\/p>","protected":false},"author":78,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[18,5,23,17],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/773"}],"collection":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/users\/78"}],"replies":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/comments?post=773"}],"version-history":[{"count":19,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/773\/revisions"}],"predecessor-version":[{"id":3419,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/773\/revisions\/3419"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/media?parent=773"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/categories?post=773"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/tags?post=773"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}