2f:T530,# Inverse iteration for quadratic Julia sets This notebook implements a `QuadraticJuliaSet` class that draws images of Julia sets using inverse iteration. ## Inverse iteration Suppose that ${tex`f:\mathbb C \to \mathbb C`} is a complex polynomial and we study the iteration of ${tex`f`}. The Julia set of ${tex`f`} can be characterized as - The boundary of the set of points that are attracted to ${tex`\infty`} or, equivalently, - The closure of the set of repelling periodic points. Either way you look at it, the Julia set is a dynamical repellor. It stands to reason, then, that the Julia set should be attractive under application of an inverse of ${tex`f`}. Inverse iteration is a strategy for generating images of Julia sets built on this observation. In the simplest and most well-known case, we have ${tex`f(z) = z^2+c`}. Of course, ${tex`f`} has two inverses: ${tex.block`f_1^{-1}(z) = \sqrt{z-c} \: \text{ and } \: f_2^{-1}(z) = -\sqrt{z-c}.`} Thus, we might try the following strategy: - Choose ${tex`z_0`} arbitrarily - Initialize a list of points ${tex`L_0 = \{z_0\}`} - Define ${tex`L_n`} recursively to contain all the inverses of all the points in ${tex`L_{n-1}`}: ${tex.block`L_n = \bigcup_{z\in L_{n-1}} \{f_1^{-1}(z),f_2^{-1}(z)\}.`} - Plot the points in ${tex`L_N`} for a suitable value of ${tex`N`}.30:T14cc,QuadraticJuliaSet = { let QuadraticJuliaSet = class QuadraticJuliaSet { constructor(c) { this.c = c; } }; QuadraticJuliaSet.prototype.draw_simple = function (opts = {}) { let { depth = 12, width = 512, height = 512, extent = [ [-2, 2], [-2, 2] ] } = opts; let c = this.c; let xmin = extent[0][0]; let xmax = extent[0][1]; let ymin = extent[1][0]; let ymax = extent[1][1]; let pts = [math.complex(0.123, -0.0456)]; for (let i = 0; i < depth; i++) { pts = pts .map(function (z) { let z1 = math.sqrt(math.add(z, math.multiply(c, -1))); let z2 = math.multiply(z1, -1); return [z1, z2]; }) .flat(); } let canvas = d3 .create("canvas") .attr("width", width) .attr("height", height); let ctx = canvas.node().getContext("2d"); ctx.fillStyle = "black"; pts = pts.map(function (z) { let ij = complex_to_canvas(z, xmin, xmax, ymin, ymax, width, height); let i = ij[0]; let j = ij[1]; ctx.fillRect(i, j, 1, 1); }); return canvas.node(); }; QuadraticJuliaSet.prototype.generate_inverse_iterate_matrix = function ( opts = {} ) { let { bailout = 5, extent = [ [-2, 2], [-2, 2] ], width = 512, height = 512 } = opts; let c = this.c; let xmin = extent[0][0]; let xmax = extent[0][1]; let ymin = extent[1][0]; let ymax = extent[1][1]; // Set up a matrix to keep track of which // points have been plotted. let plot_record = new Array(width + 1); for (let i = 0; i <= width; i++) { plot_record[i] = new Array(height + 1); } // Initialize the values in the matrix // to be zero. for (let row = 0; row <= width; row++) { for (let col = 0; col <= height; col++) { plot_record[row][col] = 0; } } // Start the iteration from an intial point. let z0 = math.complex(0.123, -0.0456); // And inverse iterate a few times to ensure // the initial point is close to the Julia set. for (let i = 0; i < 5; i++) { z0 = math.sqrt(math.add(z0, math.multiply(c, -1))); z0 = math.multiply(math.sqrt(math.add(z0, math.multiply(c, -1))), -1); } // Now, we're going to use a standard tree construction // to perform the inverse iteration. let zNode = new JuliaNode(z0); let queue = [zNode]; let cnt = 0; while (queue.length > 0) { zNode = queue.pop(); let z = zNode.z; let ij = complex_to_canvas(z, xmin, xmax, ymin, ymax, width, height); let i = ij[0]; let j = ij[1]; if (plot_record[i][j] < bailout) { let z1 = math.sqrt(math.add(z, math.multiply(c, -1))); let z2 = math.multiply(z1, -1); let left = new JuliaNode(z1); let right = new JuliaNode(z2); zNode.left = left; zNode.right = right; queue.push(left); queue.push(right); plot_record[i][j] = plot_record[i][j] + 1; } } // The output is matrix that tells us which pixels // are on (value>0) or off (value=0). return plot_record; }; QuadraticJuliaSet.prototype.draw = function (opts = {}) { let { bailout = 5, extent = [ [-2, 2], [-2, 2] ], width = 512, height = 512 } = opts; let all_opts = { bailout, extent, width, height }; let plot_record = this.generate_inverse_iterate_matrix(all_opts); let canvas = d3 .create("canvas") .attr("width", width) .attr("height", height); let ctx = canvas.node().getContext("2d"); ctx.fillStyle = "black"; let cnts = []; for (let j = 0; j <= height; j++) { for (let i = 0; i <= width; i++) { if (plot_record[i][j] > 0) { ctx.fillRect(i, j, 1, 1); } } } return canvas.node(); }; QuadraticJuliaSet.prototype.box_count = function () { return box_count( this.generate_inverse_iterate_matrix({ bailout: 5, width: 2 ** 10, height: 2 ** 10 }) ); }; return QuadraticJuliaSet; function xy_to_canvas(x, y, xmin, xmax, ymin, ymax, width, height) { return [ Math.min(Math.floor(((x - xmin) * width) / (xmax - xmin)), width - 1), Math.min(Math.floor(((ymax - y) * height) / (ymax - ymin)), height - 1) ]; } function complex_to_canvas(z, xmin, xmax, ymin, ymax, width, height) { return xy_to_canvas(z.re, z.im, xmin, xmax, ymin, ymax, width, height); } function box_count(matrix) { let w = matrix.length; let cnts = []; for (let s = 1; s <= w; s = 2 * s) { let boxes = []; let cnt = 0; for (let i = 0; s * (i + 1) < w; i++) { for (let j = 0; s * (j + 1) < w; j++) { if (check_box(i, j, s)) { cnt = cnt + 1; boxes.push([i, j]); } } } cnts.push({ s: s, cnt: cnt, boxes: boxes }); } function check_box(i, j, s) { for (let i0 = i * s; i0 < (i + 1) * s; i0++) { for (let j0 = j * s; j0 < (j + 1) * s; j0++) { if (matrix[i0][j0] > 0) { return true; } } } return false; } return cnts; } }31:[] 2d:{"title":"Inverse iteration for quadratic Julia sets","cells":[{"id":0,"value":"$2f","pinned":false,"hidden":false,"mode":"md"},{"id":393,"value":"### Using the `QuadraticJuliaSet` class\n\nThe simple inverse iterate scheme is implemented by the `draw_simple` method. We can use it to draw a Julia set as follows:","pinned":false,"hidden":false,"mode":"md"},{"id":387,"value":"new QuadraticJuliaSet(0).draw_simple()","pinned":true,"hidden":false,"mode":"ojs"},{"id":401,"value":"This illustrates the well known fact that the Julia set of ${tex`f(z)=z^2`} is the unit circle. Unfortunately, this simple technique doesn't work well for most Julia sets because the points it generates are not very well distributed throughout Julia set. The following interactive image illustrates this problem for ${tex`c=-0.62+0.41i`}. Note that the `depth` option controls the number of number of iterations for the algorithm.","pinned":false,"hidden":false,"mode":"md"},{"id":409,"value":"viewof depth = Inputs.range([1, 18], { value: 12, step: 1, label: \"Depth:\" })","pinned":false,"hidden":false,"mode":"ojs"},{"id":402,"value":"new QuadraticJuliaSet(math.complex(-0.62, 0.41)).draw_simple({ depth })","pinned":true,"hidden":false,"mode":"ojs"},{"id":446,"value":"The image above is drawn into a ${tex`512\\times512`}. At a `depth` of 18, ${tex`2^{18}`} points are generated. Now, it just so happens that\n${tex.block`2^{18} = 512 \\times 512 = ${2**18}.`}\nThus, at the highest level depth determined by the slider, we've generated enough points to fill the entire square - certainly, enough points to draw the entire Julia set at this resolution.\n\n### A modified algorithm\n\nThere's a simple modification of the inverse iteration algorithm designed to help reach those hard to get to spots. Rather than iterating until we've reached a pre-prescribed depth, we iterate from a given point only if that point's pixel has not yet been plotted a pre-prescribed number of times.\n\nHere's how the previous example looks with this modified algorithm:","pinned":false,"hidden":false,"mode":"md"},{"id":45,"value":"new QuadraticJuliaSet(math.complex(-0.62, 0.41)).draw()","pinned":true,"hidden":false,"mode":"ojs"},{"id":471,"value":"The dimension of this Julia set happens to be about 1.587 with all digits provably correct - but that's another story.","pinned":false,"hidden":false,"mode":"md"},{"id":440,"value":"## Code\n\n### The `QuadraticJuliaSet` class","pinned":false,"hidden":false,"mode":"md"},{"id":3,"value":"$30","pinned":true,"hidden":false,"mode":"ojs"},{"id":10,"value":"JuliaNode = class JuliaNode {\n constructor(z) {\n this.z = z;\n }\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":5,"value":"math = require(\"mathjs\")","pinned":true,"hidden":false,"mode":"ojs"},{"id":181,"value":"// julia_dim2 = julia_dim\n// .filter((a) => a.C < 0.2501)\n// .concat([{ C: 0.2505, Dim: undefined }])\n// .concat(julia_dim.filter((a) => a.C > 0.2505))","pinned":true,"hidden":false,"mode":"ojs"},{"id":148,"value":"// julia_dim = FileAttachment(\"julia_dim.csv\").csv({ typed: true })","pinned":true,"hidden":false,"mode":"ojs"},{"id":72,"value":"ss = require(\"simple-statistics\")","pinned":true,"hidden":false,"mode":"ojs"},{"id":531,"value":"{\n // let div = d3.create(\"div\").style(\"width\", \"600px\");\n let canvas = J.draw({\n bailout: 25,\n width: 1000,\n height: 250,\n extent: [\n [-2, 2],\n [-0.5, 0.5]\n ]\n });\n d3.select(canvas).style(\"max-width\", \"100%\");\n //div.append(() => canvas);\n return canvas;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":474,"value":"counts\n ? ss.linearRegression(counts.map((o) => [-Math.log(o.s), Math.log(o.cnt)]))\n : 0","pinned":true,"hidden":false,"mode":"ojs"},{"id":529,"value":"counts = J.box_count({\n bailout: 25,\n width: 1000,\n height: 1000\n})","pinned":true,"hidden":false,"mode":"ojs"},{"id":488,"value":"J = new QuadraticJuliaSet(c)","pinned":true,"hidden":false,"mode":"ojs"},{"id":524,"value":"// c = math.complex(-0.122561, 0.744861)\n\n// c = math.complex(0.3, 0.025)\n\nc = math.complex(-1.76, 0.005)","pinned":true,"hidden":false,"mode":"ojs"}],"files":"$W31","stdlib":"1","titling":"auto"}