31:T5ad,newton_pic = { let xmin = -2; let xmax = 2; let ymin = -1.6; let ymax = 1.6; // xmin = -0.1; // xmax = 0.3; // ymin = 0.9; // ymax = 1.2; let aspect = (ymax - ymin) / (xmax - xmin); let w = 800; let h = w * aspect; let bail = 50; let x_scale = d3 .scaleLinear() .domain([xmin, xmax]) .range([0, w]); let y_scale = d3 .scaleLinear() .domain([ymin, ymax]) .range([h, 0]); let canvas = d3 .create('canvas') .attr('width', w) .attr('height', h); let context = canvas.node().getContext("2d"); let canvasData = context.createImageData(w, h); for (let i = 0; i < w; i++) { for (let j = 0; j < h; j++) { let x = x_scale.invert(i); let y = y_scale.invert(j); let info = orbitInfo([x, y], bail); let idx = (i + j * w) * 4; let c = color(info, bail); canvasData.data[idx + 0] = c[0]; canvasData.data[idx + 1] = c[1]; canvasData.data[idx + 2] = c[2]; canvasData.data[idx + 3] = 255; } } context.putImageData(canvasData, 0, 0); context.fillStyle = 'black'; context.beginPath(); context.arc(x_scale(1), y_scale(0), 5, 0, 2 * Math.PI); context.fill(); context.beginPath(); context.arc(x_scale(-1 / 2), y_scale(Math.sqrt(3) / 2), 5, 0, 2 * Math.PI); context.fill(); context.beginPath(); context.arc(x_scale(-1 / 2), y_scale(-Math.sqrt(3) / 2), 5, 0, 2 * Math.PI); context.fill(); return canvas.node(); }32:T43c,texture = { let n = 500; let step = Math.PI / n; let theta_min = -Math.PI; let theta_max = Math.PI; let phi_min = step; let phi_max = Math.PI; // xmin = -0.1; // xmax = 0.3; // ymin = 0.9; // ymax = 1.2; let w = 2 * n; let h = n; let bail = 100; let canvas = d3.create("canvas").attr("width", w).attr("height", h); let context = canvas.node().getContext("2d"); let canvasData = context.createImageData(w, h); let j = 0; for (let theta = theta_min; theta < theta_max; theta = theta + step) { let i = 0; for (let phi = phi_min; phi < phi_max; phi = phi + step) { let d = Math.tan(phi / 2); let x = Math.cos(theta) / d; let y = Math.sin(theta) / d; let info = orbitInfo([x, y], bail); let idx = (j + i * w) * 4; let c = color(info, bail); canvasData.data[idx + 0] = c[0]; canvasData.data[idx + 1] = c[1]; canvasData.data[idx + 2] = c[2]; canvasData.data[idx + 3] = 255; i = i + 1; } j = j + 1; } context.putImageData(canvasData, 0, 0); return canvas.node(); }33:[] 2f:{"title":"Newtons's method on the Riemann sphere","cells":[{"id":0,"value":"topmatter = md`# Newtons's method on the Riemann sphere\n\nHere are the attractive basins under Newton's method applied to the polynomial ${tex`f(z) = z^3-1`} and projected onto the Riemann sphere.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":165,"value":"viewof rotation = Scrubber(rotation_values)","pinned":false,"hidden":false,"mode":"ojs"},{"id":148,"value":"riemann_sphere = {\n yield html`\n \n \n \n \n \n \n \n \n \n \n \n \n \n `;\n x3dom.reload();\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":179,"value":"md`Newton's method is a technique for finding roots of functions can lead to wonderful fractal images when applied to complex polynomials. A famous example is \n${tex.block`f(z) = z^3-1.`}\nNewton says that to find the roots of this thing, we can iterate\n${tex.block`N(z) = z - \\frac{f(z)}{f'(z)} = z - \\frac{z^3-1}{3z^2}.`}\nNow, as it turns out, the function ${tex`f`} has three roots: one real root at ${tex`z=1`} and two complex roots at\n${tex.block`z = -\\frac{1}{2} \\pm \\frac{\\sqrt{3}}{2}i.`}\nThose three points are shown as three black dots in the image below. To generate the shading, we iterate ${tex`N`} from each complex number in the rectangle and shade the corresponding pixel according to which of the three roots the process converges.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":203,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":214,"value":"md`To map onto a sphere nicely, we'll discard the rectangular grid of points that we used above in favor of a collection of points that looks something like the following (although, we'll want higher resolution, of course):\n\n![](https://i.stack.imgur.com/YElI7.png)\n\nThis grid is more nicely distributed on the sphere. For example, that the number of points inside and outside of the unit circle are the same. Using that grid of points, we can generate the following image which works nicely as a texture for our sphere:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":90,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":224,"value":"md`You can read more about this and see some more examples at [this StackExchange post](https://mathematica.stackexchange.com/questions/15047/how-to-draw-fractal-images-of-iteration-functions-on-the-riemann-sphere).`","pinned":false,"hidden":false,"mode":"ojs"},{"id":221,"value":"md`## Code`","pinned":false,"hidden":false,"mode":"ojs"},{"id":40,"value":"function color(info, bail) {\n let v = 255 * (1 - info.cnt / (bail + 1)) ** 2;\n let lim = info.lim;\n if (lim == 0) {\n return [0, 0, v];\n } else if (lim == 1) {\n return [v, 0, 0];\n } else if (lim == 2) {\n return [v, v, 0];\n } else {\n return [0, 0, 0];\n }\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":22,"value":"function orbitInfo([x0, y0], bail) {\n let cnt = 0;\n let x = x0;\n let y = y0;\n let lim = classify([x, y]);\n while (lim == 3 && cnt < bail) {\n [x, y] = N([x, y]);\n lim = classify([x, y]);\n cnt = cnt + 1;\n }\n return { lim: lim, cnt: cnt };\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":11,"value":"function classify([x, y]) {\n if ((x - 1) ** 2 + y ** 2 < 0.01) {\n return 0;\n } else if ((x + 0.5) ** 2 + (y - 0.8660254) ** 2 < 0.01) {\n return 1;\n } else if ((x + 0.5) ** 2 + (y + 0.8660254) ** 2 < 0.01) {\n return 2;\n } else {\n return 3;\n }\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":2,"value":"// The complex Newton's method iteration function expressed\n// in terms of the real and imaginary parts of z=x+iy.\nfunction N([x, y]) {\n let denominator = 3 * (x ** 2 + y ** 2) ** 2;\n let numerator_re =\n x ** 2 + 2 * x ** 5 - y ** 2 + 4 * x ** 3 * y ** 2 + 2 * x * y ** 4;\n let numerator_im =\n -2 * x * y + 2 * x ** 4 * y + 4 * x ** 2 * y ** 3 + 2 * y ** 5;\n return [numerator_re / denominator, numerator_im / denominator];\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":171,"value":"rotate = d3\n .select(riemann_sphere)\n .select('#rotation')\n .attr('rotation', `0 0 1 ${rotation}`)","pinned":true,"hidden":false,"mode":"ojs"},{"id":163,"value":"rotation_values = d3.range(0, 2 * Math.PI, Math.PI / 400)","pinned":true,"hidden":false,"mode":"ojs"},{"id":161,"value":"import { Scrubber } from \"@mbostock/scrubber\"","pinned":true,"hidden":false,"mode":"ojs"},{"id":26,"value":"d3 = require('d3@6')","pinned":true,"hidden":false,"mode":"ojs"},{"id":143,"value":"x3dom = require('x3dom').catch(() => window['x3dom'])","pinned":true,"hidden":false,"mode":"ojs"},{"id":157,"value":"html`\n \n`","pinned":true,"hidden":false,"mode":"ojs"}],"files":"$W33","stdlib":"1","titling":"auto"}