1d:{"title":"Finding Roots in the Complex Plane","cells":[{"id":3271,"value":"$1f","pinned":false,"hidden":false,"mode":"ojs"},{"id":4074,"value":"DOM.element('div') ||\n html`

\n${\n /Chrome/.test(navigator.userAgent) && /Google Inc/.test(navigator.vendor)\n ? `

🚨 Looks like you might be using Google Chrome. Due to this Chrome bug, the tall iframe that powers this page might cause lots of flickering. The solutions for now are either to temporarily disable GPU rasterization by setting the chrome://flags#enable-gpu-rasterization flag to disabled or to try the Firefox or Safari browsers instead. Simply narrowing the window also works for me.

`\n : ''\n}\n

\n\n

`","pinned":false,"hidden":false,"mode":"ojs"},{"id":4129,"value":"$20","pinned":false,"hidden":false,"mode":"ojs"},{"id":3590,"value":"md`👉Move the initial guess and observe the behavior of Newton's method.\n\n> Q: What does a root look like in this plot? (If this is unclear after a bit of experimentation, feel free to [let me know](https://twitter.com/rickyreusser)! Otherwise the article on [domain coloring](https://en.wikipedia.org/wiki/Domain_coloring) should offer assistance.)\n\n> Q: Where does it perform poorly? Where does it perform well? Why?\n\n> Q: Can we find all roots using this method?`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3386,"value":"$21","pinned":false,"hidden":false,"mode":"ojs"},{"id":3332,"value":"$22","pinned":false,"hidden":false,"mode":"ojs"},{"id":3300,"value":"newtonRaphsonPlot = createPlot({\n xrange: [-14, 14],\n yrange: [-14, 14],\n extras: {searchPoint: [6, 6]}\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":5727,"value":"html`

🌈 This notebook visualizes functions in the complex plane using domain coloring. There are two main things to know. Color represents phase (+real, -real, +imaginary, and -imaginary) while shading represents magnitude. Magnitude shading scales logarithmically so that regions with uniformly spaced contours exhibit exponential growth. When all colors meet at a point, it may indicate the presence of a zero (root). You can learn more about the general technique here and about the specific implementation here.

`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3461,"value":"md`Playing around a bit, it's clear why we're not content with Newton's method. It jumps around wildly in slowly changing regions, converges slowly in others, and doesn't consistently find the nearest root.\n\nThese are analogous to the problems Newton's method encounters in one dimension, but the extra cost and complexity of getting our bearings in two dimensions make Newton's method unacceptable in the complex plane.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":4827,"value":"md`## An aside on fractals\nJust how wildly does it jump around? Consider the function ${$$`f(z) = z^3 - 1.`} This function has three roots at ${$`1`}, ${$`-\\frac{1}{2} + \\frac{\\sqrt{3}}{2}i`}, and ${$`-\\frac{1}{2} - \\frac{\\sqrt{3}}{2}i`}. The plot below colors the complex plane by the root to which Newton iteration converges. Darker regions take more iterations to converge. The resulting pattern is called a [Newton fractal](https://en.wikipedia.org/wiki/Newton_fractal). (Tip of the hat to [Jacob Rus](https://observablehq.com/@jrus) for pointing me toward this!)\n\nThings get all complicated and lovely when we don't start right next to a root, so it perhaps suffices to say that we'd like to avoid this sort of complexity if we seek a robust root-finding algorithm.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":4913,"value":"viewof plotFractal = radio({options: [{\n label: 'root to which Newton iteration converges', value: 'root'\n}, {label: 'phase of function value, f(z)', value: 'value'}], value: 'root', description: 'color complex plane by'})","pinned":false,"hidden":false,"mode":"ojs"},{"id":4831,"value":"$23","pinned":false,"hidden":false,"mode":"ojs"},{"id":605,"value":"viewof derivatives = radio({\n options: ['Automatic differentiation', 'OES_standard_derivatives'],\n value: 'Automatic differentiation',\n description: 'Domain coloring differentiation method'\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":4836,"value":"fractalPlot = createPlot({\n xrange: [-3, 3],\n yrange: [-3, 3],\n extras: {searchPoint: [0.86, 1.6]}\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":5220,"value":"fz = ({\n glsl: `\n // Compute the value of z^3 - 1 in a glsl fragment shader. It's a vec4\n // since we use automatic differentiation to compute the derivative.\n vec4 f (vec4 z) {\n return csub(cmul(cmul(z, z), z), vec2(1, 0));\n }\n `,\n f: (function() {\n var one = [1, 0];\n var z = [0, 0, 0, 0];\n var tmp = [0, 0, 0, 0];\n // Compute the value of z^3 - 1. Unfortunately we have no choice here but\n // to reimplement the function in plain JavaScript. We once again use\n // automatic differentiation to compute the derivative, hence a little\n // bit of the complexity. You're welcome to try some different functions\n // in this cell, but until something fancier is hooked up, the Newton\n // iteration won't match the coloring unless the implementations also\n // match.\n return function f(out, a, b) {\n csetValuesD(z, a, b, 1, 0);\n csubDS(out, cmulDD(tmp, cmulDD(tmp, z, z), z), one);\n };\n })()\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":0,"value":"$24","pinned":false,"hidden":false,"mode":"ojs"},{"id":3748,"value":"md`👉 Move the contour below and observe that it counts how many roots are inside.\n\n> Q: What happens when the contour intersects a root?\n\n> Q: Can you tell how many roots there are when ${$`s_0`} is not perfectly an integer?`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3477,"value":"$25","pinned":false,"hidden":false,"mode":"ojs"},{"id":3944,"value":"cauchyPlot = createPlot({\n xrange: [-8, 8],\n yrange: [-8, 8],\n extras: {searchCircle: {z: [0, 0], r: 5}}\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":3657,"value":"md`## Problems and annoyances`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3721,"value":"md`Let's take a moment to talk about some things that can go wrong. The plot below shows the function ${$$`f(z) = \\frac{\\sqrt{z + 5 + 3i}\\sqrt{z + 5 - 3i}(z - 8 - 8i)}{(z - 3 - 7i)(z - 8 + 6i)^2}.`}\n\nFor a few different reasons, counting the zeros with ${$`s_0`} requires some care.\n\n👉Move the circle below and locate the various terms in the equation.\n\n> Q: What happens when you enclose the twice-repeated pole at ${$`8 - 6i`}?\n\n> Q: What happens when you enclose both the root and the pole in the upper-right corner?\n\n> Q: What happens when you enclose one of the square roots at ${$`-5 \\pm 3i`}? Both?`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3670,"value":"$26","pinned":false,"hidden":false,"mode":"ojs"},{"id":3959,"value":"pathologicalPlot = createPlot({\n xrange: [-16, 16],\n yrange: [-15, 15],\n extras: {searchCircle: {z: [0, 0], r: 6}}\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":3664,"value":"fPathological = ({\n glsl: `complex f (complex z) {\n return cdiv(\n cmul(\n cmul(\n csqrt(csub(z, vec2(-5, 3))),\n csqrt(csub(z, vec2(-5, -3)))\n ),\n csub(z, vec2(8, 8))\n ),\n cmul(\n csub(z, vec2(3, 7)),\n csqr(csub(z, vec2(8, -6)))\n )\n );\n }`,\n f: (function () {\n var z = ccreateD()\n var tmp1 = ccreateD()\n var tmp2 = ccreateD()\n var tmp3 = ccreateD()\n var tmp4 = ccreateD()\n return function (out, a, b) {\n csetValuesD(z, a, b, 1, 0)\n return cdivDD(out,\n cmulDD(tmp1,\n cmulDD(tmp1,\n csqrtD(tmp1, csubDS(tmp1, z, [-5, 3])),\n csqrtD(tmp2, csubDS(tmp2, z, [-5, -3]))\n ),\n csubDS(tmp2, z, [8, 8])\n ),\n cmulDD(tmp2,\n csubDS(tmp2, z, [3, 7]),\n csqrD(tmp3, csubDS(tmp4, z, [8, -6]))\n )\n );\n }\n })()\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":5256,"value":"md`The above equation has a couple things going on that can complicate our usage of the argument principle:\n- The [complex square roots](https://en.wikipedia.org/wiki/Square_root#Square_roots_of_negative_and_complex_numbers) on the left side exhibit discontinuities (\"branch cuts\") and show up in ${$`s_0`} as values of ${$`1 \\over 2`}.\n- a root and a pole together show up as ${$`s_0 = 1 - 1 = 0`}.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":4152,"value":"$27","pinned":false,"hidden":false,"mode":"ojs"},{"id":3553,"value":"md`## Recursive Subdivision\nThe ${$`s_0`} tool for counting roots suggests a simple recursive algorithm in which we just zoom in and subdivide until we isolate each individual root inside its own contour:\n\n 1. Integrate the contour and count the number of zeros by computing ${$`s_0`}.\n 2. If the region contains no zeros, stop.\n 3. Otherwise, subdivide into smaller regions.\n 4. For each subdivided region, go to step 1.\n\nThis method is simple but rather inefficient. It misses out on a lot of what complex analysis has to offer.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3865,"value":"$28","pinned":false,"hidden":false,"mode":"ojs"},{"id":4062,"value":"$29","pinned":false,"hidden":false,"mode":"ojs"},{"id":1101,"value":"function newtonsIdentities (er, ei, M, p) {\n var i, j, j2, sgn, sgn2;\n // Use Newton's Identities to construct a polynomial, the roots of\n // which match our complex analytic function:\n er[0] = 1;\n ei[0] = 0;\n er[1] = p[0];\n ei[1] = p[1];\n\n for (i = 2, sgn = -1; i <= M; i++, sgn = -sgn) {\n er[i] = p[2 * i - 2] * sgn;\n ei[i] = p[2 * i - 1] * sgn;\n for (j = 1, sgn2 = -sgn; j < i; j++, sgn2 = -sgn2) {\n j2 = 2 * (i - j - 1);\n er[i] += (er[j] * p[j2] - ei[j] * p[j2 + 1]) * sgn2;\n ei[i] += (er[j] * p[j2 + 1] + ei[j] * p[j2]) * sgn2;\n }\n er[i] /= i;\n ei[i] /= i;\n }\n\n sgn = M % 2 === 0 ? 1 : -1;\n for (i = 0; i <= M; i++, sgn = -sgn) {\n er[i] *= sgn;\n ei[i] *= sgn;\n }\n er.reverse();\n ei.reverse();\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1140,"value":"$2a","pinned":false,"hidden":false,"mode":"ojs"},{"id":3970,"value":"$2b","pinned":false,"hidden":false,"mode":"ojs"},{"id":3967,"value":"momentMethodPlot = createPlot({\n xrange: [-12, 12],\n yrange: [-10, 10],\n extras: {searchCircle: {z: [0, 0], r: 5}}\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":4509,"value":"$2c","pinned":false,"hidden":false,"mode":"ojs"},{"id":3608,"value":"$2d","pinned":false,"hidden":false,"mode":"ojs"},{"id":3953,"value":"delvesLynessPlot = createPlot({\n xrange: [-15, 15],\n yrange: [-10, 10],\n extras: {searchCircle: {z: [0, 0], r: 5}}\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":2010,"value":"viewof debugRootFinder = checkbox({options: [{label: 'Visualize search hierarchy', value: 'debug'}], value: 'debug'})","pinned":false,"hidden":false,"mode":"ojs"},{"id":2946,"value":"viewof refineRoots = checkbox({options: [{label: 'Refine roots with Newton-Raphson', value: 'refine'}], value: 'refine'})","pinned":false,"hidden":false,"mode":"ojs"},{"id":2552,"value":"viewof maxRootsPerCircle = slider({min: 1, max: 16, step: 1, value: 5, description: 'maximum number of roots per search circle'})","pinned":false,"hidden":false,"mode":"ojs"},{"id":2741,"value":"viewof maxSearchDepth = slider({min: 1, max: 16, step: 1, value: 7, description: 'maximum search depth'})","pinned":false,"hidden":false,"mode":"ojs"},{"id":2591,"value":"viewof maxIntegrationDepth = slider({min: 4, max: 20, step: 1, value: 10, description: 'adaptive integration recursion limit'})","pinned":false,"hidden":false,"mode":"ojs"},{"id":2596,"value":"viewof integrationTolerance = {\n const form = html``;\n form.i.oninput = () => form.o.value = `${(form.value = Math.pow(10, -form.i.valueAsNumber).toExponential(2))}`;\n form.i.oninput();\n return form;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":4774,"value":"$2e","pinned":false,"hidden":false,"mode":"ojs"},{"id":5311,"value":"md`## References\n- L. M. Delves, J. N. Lyness, [A Numerical Method for Locating the Zeros of an Analytic Function](https://pdfs.semanticscholar.org/1237/c0fabdd4f8c530e3b9453ce9e44e6607cedc.pdf). 1967.\n- N. I. Ioakidimis, E. G. Anastasselou, [A Modification of the Delves-Lyness Method for Locating the Zeros of Analytic Functions](https://www.sciencedirect.com/science/article/pii/0021999185901251). 1985.\n- A. Austin, P. Kravanja, L. Trefethen, [Numerical Algorithms Based on Analytic Function Values at Roots of Unity](http://people.maths.ox.ac.uk/trefethen/publication/PDF/2015_153.pdf). 2014.\n- Y. Nakatsukasa, O. Sète, L. Trefethen, [The AAA Algorithm for Rational Approximation](https://arxiv.org/pdf/1612.00337.pdf). 2017.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3157,"value":"$2f","pinned":false,"hidden":false,"mode":"ojs"},{"id":91,"value":"viewof f = slider({min: 100, max: 2e6, step: 100, value: 500e3, description: 'frequency, f, Hz'})","pinned":false,"hidden":false,"mode":"ojs"},{"id":22,"value":"viewof ν = slider({min: 0.0, max: 0.49, step: 0.01, value: 0.3, description: \"Poisson's ratio, ν\"})","pinned":false,"hidden":false,"mode":"ojs"},{"id":26,"value":"viewof θ = slider({min: 0.0, max: Math.PI, step: 0.01, value: 0.3, description: 'Viscoelastic loss angle, radians'})","pinned":false,"hidden":false,"mode":"ojs"},{"id":988,"value":"viewof parity = radio({\n options: ['Symmetric', 'Antisymmetric'],\n value: 'Symmetric',\n description: 'Lamb wave parity'\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":5286,"value":"md`Plate thickness, ${$`d`}, in meters`","pinned":false,"hidden":false,"mode":"ojs"},{"id":13,"value":"d = 10e-3 // meters","pinned":false,"hidden":false,"mode":"ojs"},{"id":5291,"value":"md`Density, ${$`\\rho`} in ${$`kg/m^3`}`","pinned":false,"hidden":false,"mode":"ojs"},{"id":9,"value":"ρ = 2700 // kg / m^3","pinned":false,"hidden":false,"mode":"ojs"},{"id":5323,"value":"md`Magnitude of Young's modulus, ${$`E`}, in ${$`Pa`}`","pinned":false,"hidden":false,"mode":"ojs"},{"id":963,"value":"Emag = 68.9e9 // Pascals","pinned":false,"hidden":false,"mode":"ojs"},{"id":5352,"value":"md`[Complex Young's modulus](https://en.wikipedia.org/wiki/Dynamic_modulus). The complex part here accounts for viscoelasticity.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1346,"value":"E = cscaleS([], [Math.cos(θ), Math.sin(θ)], Emag)","pinned":true,"hidden":false,"mode":"ojs"},{"id":5329,"value":"md`Longitudinal wave speed, ${$`c_l`}, in ${$`m/s`}. Note that of course speeds in real world must be real-valued. The imaginary part here represents spatially decaying wave modes.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1360,"value":"cl = cscaleS([], csqrtS([], E), Math.sqrt((1 - ν) / ρ / (1 + ν) / (1 - 2 * ν)))","pinned":true,"hidden":false,"mode":"ojs"},{"id":5334,"value":"md`Shear (tangential) wave speed, ${$`c_t`}, in ${$`m/s`}.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1367,"value":"ct = cscaleS([], csqrtS([], E), 1 / Math.sqrt(2 * ρ * (1 + ν)))","pinned":true,"hidden":false,"mode":"ojs"},{"id":1387,"value":"ω2 = csqrS([], cscaleS([], cinvS([], cl), 2.0 * Math.PI * f * d)) ","pinned":true,"hidden":false,"mode":"ojs"},{"id":1394,"value":"γ2 = csqrS([], cdivSS([], cl, ct))","pinned":true,"hidden":false,"mode":"ojs"},{"id":101,"value":"$30","pinned":false,"hidden":false,"mode":"ojs"},{"id":4406,"value":"md`## Plotting helpers`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3920,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":72,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":6113,"value":"import { rainbow } from '481d4902e37d39c8'","pinned":false,"hidden":false,"mode":"ojs"},{"id":6115,"value":"createRainbowTexture = regl => regl.texture(rainbow)","pinned":false,"hidden":false,"mode":"ojs"},{"id":4833,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":5489,"value":"md`These helpers in particular are highly experimental. They attempt to handle managing contexts and converting d3 scales into transformation matrices for use with WebGL, but I can't currently recommend using them!`","pinned":false,"hidden":false,"mode":"ojs"},{"id":52,"value":"import {\n reglCanvasWithOptions,\n createLayerStack,\n createViewport,\n viewportAxes,\n constrainLinearScaleAspectRatio,\n createReglViewportConfiguration,\n createReglLinearScaleConfiguration,\n createReglMap,\n persistentZoom,\n mat3ViewportFromLinearScales, mat3fromLinearScales\n} from '@rreusser/regl-tools'","pinned":false,"hidden":false,"mode":"ojs"},{"id":76,"value":"import { glslRainbow } from '@rreusser/adaptive-domain-coloring'","pinned":false,"hidden":false,"mode":"ojs"},{"id":6109,"value":"import { createPolarDomainColoringShader } from '39fe342decb60c30'","pinned":false,"hidden":false,"mode":"ojs"},{"id":4396,"value":"md`## Math helpers`","pinned":false,"hidden":false,"mode":"ojs"},{"id":4546,"value":"md`The functions imported below implement a number of complex operations with automatic differentiation. The suffix \\`D\\` refers to dual numbers with the derivative appended while \\`S\\` refers to complex scalars. Only the functions needed so far have been implemented, so the implementation is certainly incomplete.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1232,"value":"import {\n csqrD,\n ccopyD,\n csubDS,\n ccreateD,\n ccreateS,\n cmulSS,\n csqrtD,\n csinhcoshD,\n caddDD,\n cscaleD,\n cmulDD,\n csubDD,\n cfromValuesD,\n csetValuesD,\n cfromValuesS,\n cscaleS,\n csqrtS,\n cinvS,\n csqrS,\n cdivSS,\n cdivDD,\n} from '27199d72f00e07fc'","pinned":false,"hidden":false,"mode":"ojs"},{"id":3109,"value":"import {adaptiveSimpson, vectorAdaptiveSimpson} from '@rreusser/integration'","pinned":false,"hidden":false,"mode":"ojs"},{"id":1187,"value":"import { findRoots } from '@rreusser/polynomial-roots'","pinned":false,"hidden":false,"mode":"ojs"},{"id":78,"value":"import {glslComplexAutoDifferentiation} from '@rreusser/glsl-complex-auto-differentiation'","pinned":false,"hidden":false,"mode":"ojs"},{"id":50,"value":"d3 = require('d3@5')","pinned":false,"hidden":false,"mode":"ojs"},{"id":89,"value":"import {slider, radio, checkbox} from '@jashkenas/inputs'","pinned":false,"hidden":false,"mode":"ojs"},{"id":4451,"value":"$$ = tex","pinned":false,"hidden":false,"mode":"ojs"},{"id":4448,"value":"$$$ = $.block","pinned":false,"hidden":false,"mode":"ojs"},{"id":3597,"value":"html``","pinned":false,"hidden":false,"mode":"ojs"},{"id":3822,"value":"$34","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W35","stdlib":"1"}