2c:null 2d:["$","title",null,{"children":"Finding Roots in the Complex Plane | Ricky Reusser | Observable"}] 32:T8e2,md`# Finding Roots in the Complex Plane Our goal in this notebook is simple. Given a [complex function](https://en.wikipedia.org/wiki/Meromorphic_function) ${$`f(z)`}, find the *roots*, or values ${$`z_i`} at which ${$$`f(z_i) = 0.`} If ${$`f(z)`} is a polynomial (even one with complex coefficients), we're nicely set up to throw a pretty standard black box [polynomial roots function](https://github.com/scijs/durand-kerner) at it. General complex functions (like with sines or cosines, for example) are less straightforward. In fact complex analytic functions have [strict conditions on their structure](https://en.wikipedia.org/wiki/Cauchy%E2%80%93Riemann_equations) that will make this root-finding problem much more interesting and challenging than its real-valued counterpart. Let's be specific at the outset about which functions we're talking about. The techniques in this notebook won't work for just any complex function ${$`f(z)`}. If we're being precise, our complex functions must be [meromorphic](http://mathworld.wolfram.com/MeromorphicFunction.html). In very rough terms, such functions are nice at all but isolated points and include the complex analog of most of the regular functions we're used to dealing with, including addition, subtraction, multiplication, division, sines, cosines, exponentials, etc. Since this is pretty general, we won't have a difficult time coming up with them. What meromorphic functions *do not* permit is treating the real and complex parts separately, e.g. ${$`f(z) = Re(z)`}, or particularly bad ("essential") singularities. Wikipedia has [some helpful examples](https://en.wikipedia.org/wiki/Meromorphic_function#Examples) if you want to know more. Finally, I should note that I'm not a mathemetician and I haven't taken a proper course in Complex Analysis. There was a time when I needed to find roots in the complex plane to make sense of the ultrasonic waves I was measuring in the lab, but that time has passed. But here we are, and the rich and unexpected structure has made me excited to share one of my favorite problems. So if I've made mistakes, missed important things, or been overly negligent in simplifying the matter, please do [let me know](https://twitter.com/rickyreusser/status/1163332357144735744). `33:T49c,md` ## Newton's method in the complex plane We start with the first thing we always do when we want to find roots easily and efficiently: [Newton's method](https://en.wikipedia.org/wiki/Newton%27s_method). The idea of Newton's method for real-valued functions is simple: linearize the function at a point, follow that line until it hits zero, and repeat. Stated mathematically, to find the roots of function ${$`f(x)`} given an initial guess ${$`x_0`}, we repeat the iteration ${$$`x_{n+1}=x_{n}-{\frac {f(x_{n})}{f'(x_{n})}}.`} Extending the method to the complex plane is extremely simple. We replace ${$`x`} with ${$`z`} and write ${$$`z_{n+1}=z_{n}-{\frac {f(z_{n})}{f'(z_{n})}}.`} That's all there is to it! Such a simple extension doesn't work for just any two-dimensional function, but the [Cauchy-Riemann equations](https://en.wikipedia.org/wiki/Cauchy%E2%80%93Riemann_equations) which our functions obey make them very special. The plot below performs Newton's method on a function—one that will be sufficiently messy to keep this exploration interesting. (The function describes elastic wave modes in a plate. [See below](#integrationTolerance) for more details.)`34:T403,function complexNewtonRaphson (out, f, z0, tolerance, maxIterations, log) { maxIterations = maxIterations === undefined ? 50 : maxIterations tolerance = tolerance === undefined ? 1e-5 : tolerance var iteration = maxIterations var delta = 1, previousDelta = 0 var tmp = [0, 0, 0, 0] var zre = z0[0] var zim = z0[1] if (log) log.push([zre, zim]) while (delta / previousDelta > tolerance && --iteration > 0) { previousDelta = delta // Compute the function and its derivative tmp <- [f_re, f_im, f'_re, f'_im] f(tmp, zre, zim) // Divide tmp <- f / f' cdivSS(tmp, [tmp[0], tmp[1]], [tmp[2], tmp[3]]) // Update the current guess zre -= tmp[0] zim -= tmp[1] if (log) log.push([zre, zim]) // Check convergence delta = Math.hypot(tmp[0], tmp[1]) } // Return the answer if converged if (iteration > 0 || delta / previousDelta <= tolerance) { out[0] = zre out[1] = zim } else { out[0] = z0[0] out[1] = z0[1] } return out; }35:T10ef,{ const { stack, regl, xScale, yScale, currentXScale, currentYScale, container, configureRegl, viewport, searchPoint } = newtonRaphsonPlot; const svg = d3.select(stack.layers.svg); svg.selectAll('g.plot, g.axes, defs, .nrpath, g.pathPoints').remove(); const plot = svg .append('g') .attr('class', 'plot') .attr('clip-path', 'url(#clip-rect)'); var nrPath = []; var roots; const createLine = d3 .line() .x(d => currentXScale(d[0])) .y(d => currentYScale(d[1])); function computeRoots() { nrPath.length = 0; roots = [ complexNewtonRaphson( [], fLambWave.f, searchPoint.slice(), 1e-3, 100, nrPath ) ]; } const axes = svg.append('g').attr('class', 'axes'); const nrPathLine = plot .append('path') .attr('class', 'nrpath') .attr('fill', 'none') .attr('stroke', 'rgba(0, 0, 0, 0.5)') .attr('stroke-width', 2); const gPathPoints = plot.append('g').attr('class', 'pathPoints'); const gRoot = plot.append('g').attr('class', 'roots'); const handles = plot.append('g').attr('class', 'handles'); svg .append('defs') .append('clipPath') .attr('id', 'clip-rect') .append('rect') .attr('x', viewport.margin.l) .attr('y', viewport.margin.t) .attr('width', viewport.width - viewport.margin.l - viewport.margin.r) .attr('height', viewport.height - viewport.margin.t - viewport.margin.b); function updateRoots() { gRoot.selectAll('g.root').remove(); const enter = gRoot .selectAll('g.root') .data(roots) .enter() .append('g') .attr('class', 'root'); enter .append('circle') .attr('r', 6) .attr('stroke-width', 0) .attr('stroke', 'transparent') .attr('fill', 'rgba(255, 255, 255, 1)') .attr('cx', d => currentXScale(d[0])) .attr('cy', d => currentYScale(d[1])); enter .append('circle') .attr('r', 4) .attr('stroke-width', 0) .attr('stroke', 'transparent') .attr('fill', 'rgba(0, 0, 0, 1)'); enter .selectAll('circle') .attr('cx', d => currentXScale(d[0])) .attr('cy', d => currentYScale(d[1])); handles .selectAll('circle') .attr('cx', d => currentXScale(searchPoint[0])) .attr('cy', d => currentYScale(searchPoint[1])); gPathPoints.selectAll('circle').remove(); gPathPoints .selectAll('circle') .data(nrPath) .enter() .append('circle') .attr('fill', 'rgba(255, 255, 255, 1)') .attr('r', 2.5) .attr('stroke-width', 1) .attr('stroke', 'rgba(0, 0, 0, 1)') .attr('cx', d => currentXScale(d[0])) .attr('cy', d => currentYScale(d[1])); nrPathLine.datum(nrPath).attr("d", createLine); } function updateAxes() { axes.call(viewportAxes(viewport, currentXScale, currentYScale)); } const drawFunction = createDrawFunction(regl, fLambWave); function updatePlot() { updateRoots(); regl.poll(); regl.clear({ color: [1, 1, 1, 1] }); configureRegl.viewport(viewport, () => { configureRegl.scale(currentXScale, currentYScale, () => { configureRegl.map(() => { drawFunction({ ω2, γ2 }); }); }); }); } handles .selectAll('g.handles') .data(['z0', 'r0']) .enter() .append('g') .attr('class', 'handle') .call(el => { el.append('circle') .attr('r', 6) .attr('stroke-width', 2) .attr('stroke', 'rgba(255, 255, 255, 0.9)') .attr('fill', '#35e'); el.append('circle') .attr('cursor', 'move') .attr('r', 40) .attr('stroke-width', 0) .attr('stroke', 'transparent') .attr('fill', 'transparent') .call( d3.drag().on('drag', d => { searchPoint[0] = currentXScale.invert(d3.event.x); searchPoint[1] = currentYScale.invert(d3.event.y); computeRoots(); updateRoots(); }) ); }); svg.call( persistentZoom(currentXScale, currentYScale, xScale, yScale) .scaleExtent([1, 1]) .on('zoom.axes', updateAxes) .on('zoom.plot', updatePlot) ); computeRoots(); updateAxes(); updatePlot(); updateRoots(); yield container; }36:T1092,{ const {stack, regl, xScale, yScale, currentXScale, currentYScale, container, configureRegl, viewport, searchPoint} = fractalPlot; const func = fz const drawFunction = (plotFractal === 'root' ? createFractalDrawFunction : createDrawFunction)(regl, func) const svg = d3.select(stack.layers.svg) svg.selectAll('g.plot, g.axes, defs, .nrpath, g.pathPoints').remove(); const plot = svg.append('g') .attr('class', 'plot') .attr('clip-path', 'url(#clip-rect)'); var nrPath = [] var roots const createLine = d3.line() .x(d => currentXScale(d[0])) .y(d => currentYScale(d[1])); function computeRoots () { nrPath.length = 0 roots = [complexNewtonRaphson([], func.f, searchPoint.slice(), 1e-3, 100, nrPath)] } const axes = svg.append('g').attr('class', 'axes'); const nrPathLine = plot.append('path') .attr('class', 'nrpath') .attr('fill', 'none') .attr('stroke', 'rgba(0, 0, 0, 0.5)') .attr('stroke-width', 2) const gPathPoints = plot.append('g').attr('class', 'pathPoints'); const gRoot = plot.append('g').attr('class', 'roots'); const handles = plot.append('g').attr('class', 'handles') svg.append('defs') .append('clipPath').attr('id', 'clip-rect') .append('rect') .attr('x', viewport.margin.l) .attr('y', viewport.margin.t) .attr('width', viewport.width - viewport.margin.l - viewport.margin.r) .attr('height', viewport.height - viewport.margin.t - viewport.margin.b); function updateRoots () { gRoot.selectAll('g.root').remove() const enter = gRoot.selectAll('g.root') .data(roots) .enter().append('g').attr('class', 'root') enter.append('circle') .attr('r', 6) .attr('stroke-width', 0) .attr('stroke', 'transparent') .attr('fill', 'rgba(255, 255, 255, 1)') .attr('cx', d => currentXScale(d[0])) .attr('cy', d => currentYScale(d[1])) enter.append('circle') .attr('r', 4) .attr('stroke-width', 0) .attr('stroke', 'transparent') .attr('fill', 'rgba(0, 0, 0, 1)') enter.selectAll('circle') .attr('cx', d => currentXScale(d[0])) .attr('cy', d => currentYScale(d[1])) handles.selectAll('circle') .attr('cx', d => currentXScale(searchPoint[0])) .attr('cy', d => currentYScale(searchPoint[1])) gPathPoints.selectAll('circle').remove() gPathPoints.selectAll('circle').data(nrPath) .enter().append('circle') .attr('fill', 'rgba(255, 255, 255, 1)') .attr('r', 2.5) .attr('stroke-width', 1) .attr('stroke', 'rgba(0, 0, 0, 1)') .attr('cx', d => currentXScale(d[0])) .attr('cy', d => currentYScale(d[1])) nrPathLine.datum(nrPath).attr("d", createLine); } function updateAxes () { axes.call(viewportAxes(viewport, currentXScale, currentYScale)) } function updatePlot () { updateRoots() regl.poll() regl.clear({color: [1, 1, 1, 1]}) configureRegl.viewport(viewport, () => { configureRegl.scale(currentXScale, currentYScale, () => { configureRegl.map(() => { drawFunction({ω2, γ2}) }) }) }) } handles.selectAll('g.handles') .data(['z0', 'r0']) .enter().append('g') .attr('class', 'handle') .call(el => { el.append('circle') .attr('r', 6) .attr('stroke-width', 2) .attr('stroke', 'rgba(255, 255, 255, 0.9)') .attr('fill', '#35e') el.append('circle') .attr('cursor', 'move') .attr('r', 40) .attr('stroke-width', 0) .attr('stroke', 'transparent') .attr('fill', 'transparent') .call(d3.drag().on('drag', d => { searchPoint[0] = currentXScale.invert(d3.event.x) searchPoint[1] = currentYScale.invert(d3.event.y) computeRoots() updateRoots() })); }) svg.call(persistentZoom(currentXScale, currentYScale, xScale, yScale) .scaleExtent([0.01, 100000]) .on('zoom.axes', updateAxes) .on('zoom.plot', updatePlot)) computeRoots(); updateAxes(); updatePlot(); updateRoots() yield container }37:Tc71,md`## Cauchy's Argument Principle We seek a better method. The most memorable theorem in what I did learn of Complex Analysis is [Cauchy's Argument Principle](https://en.wikipedia.org/wiki/Argument_principle). It states that upon integrating a complex function ${$`f(z)`} around a closed contour ${$`C`}, we have${$$`\displaystyle s_0 \equiv \frac{1}{2 \pi i}\oint _{C}{f'(z) \over f(z)}\,dz=Z-P,`} where ${$`Z`} and ${$`P`} are the number of _zeros_ and _poles_ fully enclosed by the contour, respectively. We label this quantity ${$`s_0`} for reasons that will become apparent. Let's tease this apart slowly. *Zeros* (which we somewhat interchangeably call *roots*) are points where the function has the value ${$`0 + 0i`}. *Poles* are points where we divide by zero and the magnitude diverges to infinity. For the sort of function we're looking at, zeros and poles may only occur at isolated points rather than over whole regions. Working through the terms of ${$`s_0`} one by one, the first contains the number ${$`2\pi i`}, which is just a constant [complex number](https://en.wikipedia.org/wiki/Complex_number) by which we divide. Next, the quantity we integrate, ${$`f'(z) \over f(z)`}, is the derivative of the complex function at point ${$`z`} divided by the value. That means we need the derivative, which can be a little tricky. [People try hard](https://www.sciencedirect.com/science/article/pii/0021999185901251) to avoid needing the derivative, but I've used [automatic differentiation](https://en.wikipedia.org/wiki/Automatic_differentiation) in this notebook (both [in JavaScript](https://observablehq.com/@rreusser/complex-auto-differentiation) and [in GLSL](https://observablehq.com/@rreusser/glsl-complex-auto-differentiation)) to compute it automatically along with the function value. Finally, the symbol ${$`\oint_C \cdots dz`} represents an [integral in the complex plane](https://en.wikipedia.org/wiki/Contour_integration). Just like integration of real numbers, it represents the sum of lots of little pieces, except instead of two limits of integration along the real number line (a lower and upper limit, commonly the ${$`a`} and ${$`b`} in ${$`\int_a^b f(x) dx`}), we choose a path, or [*contour*](https://en.wikipedia.org/wiki/Contour_integration#Contours), in the two-dimensional complex plane. Each small step along the path is represented by the complex differential ${$`dz`}, and the circle in the integral symbol ${$`\oint`} means the path is closed. In short, this somewhat ugly-looking integral is just a value we can compute relatively easily by stepping around a path and adding up the little pieces of the derivative divided by the value. Upon doing so, we find that we've counted the difference between the number of zeros and the number of poles! It's not as unlikely as one may think to have encountered this in real science and engineering applications—or perhaps in 3Blue1Brown's excellent video "[Solving 2D equations using color, a story of winding numbers and composition](https://www.youtube.com/watch?v=b7FxPsqfkOY)." In fact, if the above equation doesn't make much sense, this would be a great time to go watch that video!` 38:T1180,{ const {stack, regl, xScale, yScale, currentXScale, currentYScale, container, configureRegl, viewport, searchCircle} = cauchyPlot; const svg = d3.select(stack.layers.svg) svg.selectAll('g.plot, g.axes, defs').remove(); const plot = svg.append('g') .attr('class', 'plot') .attr('clip-path', 'url(#clip-rect)'); var searchOutput = [] var roots function computeRoots () { roots = findComplexRoots(fLambWave.f, searchCircle.z, searchCircle.r, { s0IntegrationTolerance: 1e-8, s0IntegrationDepth: 8, integrationTolerance: 1e-10, maxIntegrationDepth: 1, searchOutput, maxRootsPerCircle: 1, maxSearchDepth: 1 }) } let axes = svg.append('g').attr('class', 'axes'); const gRoot = plot.append('g').attr('class', 'roots'); const gSearchCircles = plot.append('g').attr('class', 'searchCircles'); let handles = plot.append('g').attr('class', 'handles') svg.append('defs') .append('clipPath').attr('id', 'clip-rect') .append('rect') .attr('x', viewport.margin.l) .attr('y', viewport.margin.t) .attr('width', viewport.width - viewport.margin.l - viewport.margin.r) .attr('height', viewport.height - viewport.margin.t - viewport.margin.b); let searchCircles function updateRoots () { let searchCircleSelection = gSearchCircles.selectAll('g.searchCircle') .data(searchOutput) .join(enter => { var appended = enter.append('g') .attr('class', 'searchCircle') appended.append('circle') .attr('stroke-width', 1) .attr('stroke', 'rgba(0, 0, 0, 0.8)') .attr('fill', 'none'); appended.append('text') .attr("font-style", "italic") .attr("font-size", "20px") .attr('dy', 47) .attr('alignment-baseline', 'text-top') .attr('text-anchor', 'middle') .attr('stroke', 'rgba(255, 255, 255, 0.7)') .attr('stroke-width', 3) .attr('paint-order', 'stroke') .attr('stroke-linecap', 'butt') .attr('stroke-linejoin', 'miter') .attr('stroke-miterlimit', 2); }) searchCircleSelection.exit().remove() gSearchCircles.selectAll('g.searchCircle') .select('circle') .attr('cx', d => currentXScale(d.z[0])) .attr('cy', d => currentYScale(d.z[1])) .attr('r', d => currentXScale(d.r) - currentXScale(0)) .attr('stroke', d => 'rgba(0, 0, 0, 1.0)') gSearchCircles.selectAll('g.searchCircle').select('text') .attr('x', d => currentXScale(d.z[0])) .attr('y', d => currentYScale(d.z[1])) .text(d => d.text) handles.selectAll('circle') .attr('cx', d => currentXScale(searchCircle.z[0] + (d === 'z0' ? 0 : searchCircle.r))) .attr('cy', d => currentYScale(searchCircle.z[1])) } function updateAxes () { axes.call(viewportAxes(viewport, currentXScale, currentYScale)) } const drawFunction = createDrawFunction(regl, fLambWave) function updatePlot () { updateRoots() regl.poll() regl.clear({color: [1, 1, 1, 1]}) configureRegl.viewport(viewport, () => { configureRegl.scale(currentXScale, currentYScale, () => { configureRegl.map(() => { drawFunction({ω2, γ2}) }) }) }) } handles.selectAll('g.handles') .data(['z0', 'r0']) .enter().append('g') .attr('class', 'handle') .call(el => { el.append('circle') .attr('r', 6) .attr('stroke-width', 2) .attr('stroke', 'rgba(255, 255, 255, 0.9)') .attr('fill', '#35e') el.append('circle') .attr('cursor', 'move') .attr('r', 40) .attr('stroke-width', 0) .attr('stroke', 'transparent') .attr('fill', 'transparent') .call(d3.drag().on('drag', d => { if (d === 'r0') { searchCircle.r = currentXScale.invert(d3.event.x) - searchCircle.z[0]; } else { searchCircle.z[0] = currentXScale.invert(d3.event.x) searchCircle.z[1] = currentYScale.invert(d3.event.y) } computeRoots() updatePlot() })); }) svg.call(persistentZoom(currentXScale, currentYScale, xScale, yScale) .scaleExtent([1, 1]) .on('zoom.axes', updateAxes) .on('zoom.plot', updatePlot)) computeRoots() updateAxes(); updatePlot(); updateRoots() yield cauchyPlot.container }39:T11d8,{ const {stack, regl, xScale, yScale, currentXScale, currentYScale, container, configureRegl, viewport, searchCircle} = pathologicalPlot; const svg = d3.select(stack.layers.svg) svg.selectAll('g.plot, g.axes, defs').remove(); const plot = svg.append('g') .attr('class', 'plot') .attr('clip-path', 'url(#clip-rect)'); var searchOutput = [] var roots function computeRoots () { roots = findComplexRoots(fPathological.f, searchCircle.z, searchCircle.r, { s0IntegrationTolerance: 1e-9, s0IntegrationDepth: 11, integrationTolerance: 1e-10, maxIntegrationDepth: 1, searchOutput, maxRootsPerCircle: 1, maxSearchDepth: 1 }) } let axes = svg.append('g').attr('class', 'axes'); const gRoot = plot.append('g').attr('class', 'roots'); const gSearchCircles = plot.append('g').attr('class', 'searchCircles'); let handles = plot.append('g').attr('class', 'handles') svg.append('defs') .append('clipPath').attr('id', 'clip-rect') .append('rect') .attr('x', viewport.margin.l) .attr('y', viewport.margin.t) .attr('width', viewport.width - viewport.margin.l - viewport.margin.r) .attr('height', viewport.height - viewport.margin.t - viewport.margin.b); let searchCircles function updateRoots () { let searchCircleSelection = gSearchCircles.selectAll('g.searchCircle') .data(searchOutput) .join(enter => { var appended = enter.append('g') .attr('class', 'searchCircle') appended.append('circle') .attr('stroke-width', 1) .attr('stroke', 'rgba(0, 0, 0, 0.8)') .attr('fill', 'none'); appended.append('text') .attr("font-style", "italic") .attr("font-size", "20px") .attr('dy', 47) .attr('alignment-baseline', 'text-top') .attr('text-anchor', 'middle') .attr('stroke', 'rgba(255, 255, 255, 0.7)') .attr('stroke-width', 3) .attr('paint-order', 'stroke fill') .attr('stroke-linecap', 'butt') .attr('stroke-linejoin', 'miter') .attr('stroke-miterlimit', 2); }) searchCircleSelection.exit().remove() gSearchCircles.selectAll('g.searchCircle') .select('circle') .attr('cx', d => currentXScale(d.z[0])) .attr('cy', d => currentYScale(d.z[1])) .attr('r', d => currentXScale(d.r) - currentXScale(0)) .attr('stroke', d => 'rgba(0, 0, 0, 1.0)') gSearchCircles.selectAll('g.searchCircle').select('text') .attr('x', d => currentXScale(d.z[0])) .attr('y', d => currentYScale(d.z[1])) .text(d => d.text) handles.selectAll('circle') .attr('cx', d => currentXScale(searchCircle.z[0] + (d === 'z0' ? 0 : searchCircle.r))) .attr('cy', d => currentYScale(searchCircle.z[1])) } function updateAxes () { axes.call(viewportAxes(viewport, currentXScale, currentYScale)) } const drawFunction = createDrawFunction(regl, fPathological) function updatePlot () { updateRoots() } function updateFunction () { regl.poll() regl.clear({color: [1, 1, 1, 1]}) configureRegl.viewport(viewport, () => { configureRegl.scale(currentXScale, currentYScale, () => { configureRegl.map(() => { drawFunction({ω2, γ2}) }) }) }) } handles.selectAll('g.handles') .data(['z0', 'r0']) .enter().append('g') .attr('class', 'handle') .call(el => {el.append('circle') .attr('r', 7) .attr('stroke-width', 0) .attr('stroke', 'transparent') .attr('fill', 'rgba(255, 255, 255, 0.9)') el.append('circle') .attr('cursor', 'move') .attr('r', 5) .attr('stroke-width', 30) .attr('stroke', 'transparent') .attr('fill', '#35e') .call(d3.drag().on('drag', d => { if (d === 'r0') { searchCircle.r = currentXScale.invert(d3.event.x) - searchCircle.z[0]; } else { searchCircle.z[0] = currentXScale.invert(d3.event.x) searchCircle.z[1] = currentYScale.invert(d3.event.y) } computeRoots() updatePlot() })); }) svg.call(persistentZoom(currentXScale, currentYScale, xScale, yScale) .scaleExtent([1, 1]) .on('zoom.axes', updateAxes) .on('zoom.plot', updatePlot) .on('zoom.fn', updateFunction)) computeRoots() updateAxes(); updatePlot(); updateRoots(); updateFunction() yield container }3a:T63a,md`There are a few ways these difficulties can play out. If ${$`s_0`} is near an integer, we can hope it correctly counts the number of enclosed zeros ${$`Z`}—but we can't be certain. When we enclose both terms with a square root, ${$`s_0 = \frac{1}{2} + \frac{1}{2} = 1`}, and there's not immediately a way to tell whether it's two square roots or a single zero (or two zeros and a pole, for that matter). It's my (limited) experience that the best way around many of these difficulties is to try avoiding them to begin with. For example, consider solutions of the equation ${$$`\tan z = 1.`} We can cast solutions as the root-finding problem ${$$`f(z) = \tan z - 1 = 0,`} but this equation has poles (division by zero) that make life with our somewhat elementary tools more difficult. Instead, we may write ${$$`\frac{\sin z}{\cos z} - 1 = 0`} and multiply through by ${$`\cos z`} to obtain ${$$`f(z) = \sin z - \cos z = 0.`} Unlike ${$`\tan z`} which has poles that may confuse our interpretation of ${$`s_0,`} this equivalent equation has only roots, so that rather than figuring out how to interpret poles, we've avoided them to begin with. Branch cuts (discontinuities resulting from things like square roots) that result in non-integer values of ${$`s_0`} are a bit harder to deal with but can sometimes be worked around by rearrangement or by avoiding contours large enough to include multiple branch cuts. Of course these workarounds don't work in every situation, but they're sufficiently useful that for the rest of this notebook we'll neglect some of these thornier problems.`3b:T80c,md`## The Method of Delves and Lyness Complex analysis has been around for a while, so it's not surprising the work on numerical methods goes pretty far back. This whole notebook does little more than recite the development from Delves' and Lyness' 1967 paper, [A Numerical Method for Locating the Zeros of an Analytic Function](https://pdfs.semanticscholar.org/1237/c0fabdd4f8c530e3b9453ce9e44e6607cedc.pdf). In their paper, Delves and Lyness note as a generalization of Cauchy's Argument Principle the fact that ${$$`\displaystyle s_N = \frac{1}{2 \pi i}\oint_{C}z^N{f'(z) \over f(z)}\,dz=\sum \limits_{i = 1}^{Z} (z_i)^N`} where ${$`z_i\;(i = 1,\,2,\,\cdots,\,Z)`} are the zeros inside the region bounded by ${$`C`}. Another ugly equation. First let's consider the particular case of ${$`N = 0`}. Since ${$`z^0 = 1`}, we have ${$$` \begin{aligned} s_0 =& \frac{1}{2 \pi i}\oint _{C}{f'(z) \over f(z)}\,dz \\ =& \sum \limits_{i = 1}^{Z} (z_i)^0 = \sum \limits_{i = 1}^{Z} 1 = Z, \end{aligned} `} thereby reducing to our previous formula for counting the zeros (allowing for the time being that we may neglect poles). For other ${$`N`}, it's at the very least a value we may compute directly by way of numerical integration. So basically ${$`s_N`} are like ${$`s_0`} but with powers of ${$`z`} stuck inside. The important thing here is that the formula doesn't just allow computation of one value. We can plug in various ${$`N`} and—with a bit of straightforward numerical integration—compute as many ${$`s_N`} as we'd like. Therefore considering ${$`s_N`} as good as computed numbers, we relate them to the zeros via the formula and write ${$$`\displaystyle \begin{aligned} s_1 &= \sum \limits_{i = 1}^{Z} z_i^1 = z_1 + z_2 + \cdots + z_Z \\ s_2 &= \sum \limits_{i = 1}^{Z} z_i^2 = z_1^2 + z_2^2 + \cdots + z_Z^2 \\ & \vdots \\ s_Z &= \sum \limits_{i = 1}^{Z} z_i^Z = z_1^Z + z_2^Z + \cdots + z_Z^Z \end{aligned} `} This looks an awful lot like a system of ${$`Z`} equations with ${$`Z`} unknowns which usually indicates a solved problem.`3c:Tb77,md`## Symmetric polynomials and Newton's identities Despite the apparent simplicity of our ${$`Z`} equations with ${$`Z`} unknowns, they're not linear so that we require more than our normal linear algebra tools to solve them. Equations of the form ${$$`s_k(z_1, z_2, \ldots,z_N) = \sum_{i = 1}^Nz_i^k`} are called [power sum symmetric polynomials](https://en.wikipedia.org/wiki/Power_sum_symmetric_polynomial). They're called symmetric because swapping any of the ${$`z_i`} leaves ${$`s`} unchanged. In other words, we don't care which root is which; we only care how they behave in aggregate. Symmetric polynomials are fascinating objects though are certainly not a topic I'm qualified to write about. More than anything, I've been surprised to encounter them in a few different entirely applied contexts, my favorite of which is Mikola Lysenko's delightful article, "[Comparing Sequences Without Sorting](https://0fps.net/2013/01/24/comparing-sequences-without-sorting/)" in which he makes use of the fact that symmetric polynomials have all the right properties—in particular, not caring which is which—to compare two sequences of numbers without first sorting them. To make a long story short, [Newton's identities](https://en.wikipedia.org/wiki/Newton%27s_identities) give us the inversion we seek in the form of a recursive definition. In terms of the elementary symmetric polynomials ${$`e_k`} (about which details we need not worry here), we write ${$$` \begin{aligned} e_{0} &=1, \\ -e_{1}&=-s_{1}, \\ e_{2} &={\frac {1}{2}}(e_{1}s_{1}-s_{2}), \\ -e_{3}&=-{\frac {1}{3}}(e_{2}s_{1}-e_{1}s_{2}+s_{3}), \\ e_{4} &={\frac {1}{4}}(e_{3}s_{1}-e_{2}s_{2}+e_{1}s_{3}-s_{4}), \\ & \vdots \end{aligned} `} together with the polynomial ${$$` \prod _{i=1}^{Z}\left(z-z_{i}\right)=\sum _{k=0}^{Z}(-1)^{k}e_{k}z^{Z-k}. `} The lefthand side, by virtue of its definition, is a polynomial containing the unknown zeros we're trying to find. The righthand side is a polynomial we may compute by plugging our known ${$`s_N`} into the inversion Newton's identities gives us. From there, any [lovely polynomial root-finder](https://github.com/scijs/durand-kerner) will complete the job and yield the roots of the complex function! Restating this all in a much simpler manner, *the method of Delves and Lyness converts the problem of finding roots in the *interior* of a closed region to a problem of analyzing values on the *boundary* of that region*. What's most amazing to me is that the computed polynomial contains the roots of ${$`f(z)`} even though the function itself doesn't have to be any sort of simple polynomial. The plot below implements this method. For the circular contour, ${$`s_0 = Z`} values of ${$`s_1, s_2, \ldots, s_Z`} are computed and used to construct a polynomial whose roots coincide with the roots inside the contour. Observe that if you ensnare enough roots, it gets pretty unstable or fails entirely. `3d:T1b9d,findComplexRoots = { var out = new Array(2); var params = {}; function cmod (a, b) { var c var aa = Math.abs(a) var ab = Math.abs(b) if (aa === 0 && ab === 0) return 0 if (aa >= ab) { c = b / a return aa * Math.sqrt(1 + c * c) } c = a / b return ab * Math.sqrt(1 + c * c) } function cdiv (out, a, b, c, d ) { var e, f if (Math.abs(c) >= Math.abs(d)) { e = d / c f = 1 / (c + d * e) out[0] = (a + b * e) * f out[1] = (b - a * e) * f } else { e = c / d f = 1 / (c * e + d) out[0] = (a * e + b) * f out[1] = (b * e - a) * f } return out } var dzrsub = [ 0, 0.81, 0.5727564927611035, 0, -0.5727564927611035, -0.81, -0.5727564927611035, 0, 0.5727564927611035 ]; var dzisub = [ 0, 0, 0.5727564927611035, 0.81, 0.5727564927611035, 0, -0.5727564927611035, -0.81, -0.5727564927611035 ]; var rsub = [ 0.5, 0.4166666666666667, 0.4166666666666667, 0.4166666666666667, 0.4166666666666667, 0.4166666666666667, 0.4166666666666667, 0.4166666666666667, 0.4166666666666667 ]; function s0Integrand (theta) { var c, s, a, b, rtmp, itmp; c = Math.cos(theta); s = Math.sin(theta); a = params.z0r + c * params.r; b = params.z0i + s * params.r; // tmp <- f(z) params.f(out, a, b); rtmp = out[0]; itmp = out[1]; // out <- f'(z) out[0] = out[2]; out[1] = out[3]; // out /= f(z) cdiv(out, out[0], out[1], rtmp, itmp); // re(tmp * dz / (2 pi i)) return out[0] * c - out[1] * s; } function smIntegrand (out, theta) { var c, s, a, b, rtmp, itmp; c = Math.cos(theta); s = Math.sin(theta); // out <- f(z) params.f(out, params.z0r + c * params.r, params.z0i + s * params.r); // out <- f'(z) / f(z) cdiv(out, out[2], out[3], out[0], out[1]); // out *= dz (/ (2 pi) but leave this part until outside the loop) rtmp = out[0]; out[0] = out[0] * c - out[1] * s; out[1] = rtmp * s + out[1] * c; // out *= z rtmp = out[0]; out[0] = (out[0] * c - out[1] * s); out[1] = (rtmp * s + out[1] * c); // Recursively multiply by z to compute successive moments z^N: // out[i] = out[i - 1] * z for (var i = 2; i < 2 * params.M; i += 2) { out[i] = out[i - 2] * c - out[i - 1] * s; out[i + 1] = out[i - 2] * s + out[i - 1] * c; } return out } function locateZeros (qr, qi, f, zr, zi, r, z0r, z0i, r0, tol, maxIntegrationDepth, s0IntegrationTolerance, s0IntegrationDepth, curDepth, maxRootsPerCircle, searchCircleOutput) { var p, s0, er, ei, i, j, qri, qii, qmod, qC, M, rr; var qNr, qNi, tmp, qNew, success; // Exit if we've reached the recursion limit: if (curDepth <= 0) { return false; } // Exit early if this circle is completely outside the search region if (cmod(zr - z0r, zi - z0i) - r > r0) { return true } params.z0r = zr; params.z0i = zi; params.r = r; params.f = f; s0 = adaptiveSimpson(s0Integrand, 0, Math.PI * 2, s0IntegrationTolerance, s0IntegrationDepth) * r * (0.5 / Math.PI); M = Math.round(s0); qC = 0; for (i = qr.length - 1; i >= 0; i--) { // Count this zero if it's inside the circle: if (cmod(qr[i] - zr, qi[i] - zi) < r) qC++; } if (M - qC <= 0.9 || s0 < 0.9) { } else if (Math.abs(s0 - M) > 1e-1 || M - qC > maxRootsPerCircle) { // If this happens, then we had a pretty inaccurate integration: for (i = 0, success = true; i < 9; i++) { rr = r * rsub[i]; success = locateZeros(qr, qi, f, zr + r * dzrsub[i], zi + r * dzisub[i], rr, z0r, z0i, r0, tol, maxIntegrationDepth, s0IntegrationTolerance, s0IntegrationDepth, curDepth - 1, maxRootsPerCircle, searchCircleOutput ) && success; } } else { // Compute the zeros: params.M = M - qC; var p = vectorAdaptiveSimpson([], smIntegrand, 0, Math.PI * 2, tol, maxIntegrationDepth, 2 * params.M); // post-multiply the integral by this constant factor var f = r * 0.5 / Math.PI; for (var i = 2 * params.M - 1; i >= 0; i--) { p[i] *= f; } // Deflate the result by known zeros inside the current circle: for (i = qr.length - 1; i >= 0; i--) { // Scale to the unit circle in which we're working qri = (qr[i] - zr) / r; qii = (qi[i] - zi) / r; if (cmod(qri, qii) < 1) { // Initialize the powers of qN var qNr = qri; var qNi = qii; for (j = 0, qNr = 1, qNi = 0; j < p.length; j += 2) { tmp = qNr; qNr = tmp * qri - qNi * qii; qNi = tmp * qii + qNi * qri; // Subtract off the power of the zero: p[j] -= qNr; p[j + 1] -= qNi; } } } // Apply Newton's Identities, constructing a polynomial sharing the roots // of the analytic function: er = []; ei = []; newtonsIdentities(er, ei, params.M, p); // Compute the roots of the constructed polynomial: qNew = findRoots(er, ei, 100 * params.M * params.M, tol); // Append and filter to zeros inside the original contour: for (i = qNew[0].length - 1; i >= 0; i--) { var qrRescaled = zr + qNew[0][i] * r; var qiRescaled = zi + qNew[1][i] * r; //if (cmod(qrRescaled - z0r, qiRescaled - z0i) < r0) { qr.push(qrRescaled); qi.push(qiRescaled); //} } } if (searchCircleOutput) { searchCircleOutput.push({z: [zr, zi], r: r, text: 's₀ = ' + s0.toFixed(5).toString().replace(/(\.[0-9]*[1-9]|\.0)0*$/,'$1').replace(/\.0$/,'').replace(/^-0$/, '0')}) } return [qr, qi]; } return function delvesLyness (f, z0, r, opts) { opts = opts || {} var s0IntegrationDepth = opts.s0IntegrationDepth === undefined ? 13 : opts.s0IntegrationDepth var s0IntegrationTolerance = opts.s0IntegrationTolerance === undefined ? 1e-6 : opts.s0IntegrationTolerance var maxIntegrationDepth = opts.maxIntegrationDepth === undefined ? 12 : opts.maxIntegrationDepth var integrationTolerance = opts.integrationTolerance === undefined ? 1e-8 : opts.integrationTolerance var maxSearchDepth = opts.maxSearchDepth === undefined ? 4 : opts.maxSearchDepth var maxRootsPerCircle = opts.maxRootsPerCircle === undefined ? 4 : opts.maxRootsPerCircle if (opts.searchOutput) opts.searchOutput.length = 0; var foundRoots = locateZeros([], [], f, z0[0], z0[1], r, z0[0], z0[1], r, integrationTolerance, maxIntegrationDepth, s0IntegrationTolerance, s0IntegrationDepth, maxSearchDepth, maxRootsPerCircle, opts.searchOutput) var output = []; for (var i = 0; i < foundRoots[0].length; i++) { if (cmod(foundRoots[0][i] - z0[0], foundRoots[1][i] - z0[1]) <= r) { output.push([foundRoots[0][i], foundRoots[1][i]]) } } return output; } }3e:T14af,{ const {stack, regl, xScale, yScale, currentXScale, currentYScale, container, configureRegl, viewport, searchCircle} = momentMethodPlot; const svg = d3.select(stack.layers.svg) svg.selectAll('g.plot, g.axes, defs').remove(); const plot = svg.append('g') .attr('class', 'plot') .attr('clip-path', 'url(#clip-rect)'); var searchOutput = [] var roots function computeRoots () { roots = findComplexRoots(fLambWave.f, searchCircle.z, searchCircle.r, { integrationTolerance: 1e-11, maxIntegrationDepth: 11, searchOutput, maxRootsPerCircle: Infinity, maxSearchDepth: 1, }) } let axes = svg.append('g').attr('class', 'axes'); const gRoot = plot.append('g').attr('class', 'roots'); const gSearchCircles = plot.append('g').attr('class', 'searchCircles'); let handles = plot.append('g').attr('class', 'handles') svg.append('defs') .append('clipPath').attr('id', 'clip-rect') .append('rect') .attr('x', viewport.margin.l) .attr('y', viewport.margin.t) .attr('width', viewport.width - viewport.margin.l - viewport.margin.r) .attr('height', viewport.height - viewport.margin.t - viewport.margin.b); let searchCircles function updateRoots () { gRoot.selectAll('g.root').remove() const enter = gRoot.selectAll('g.root') .data(roots) .enter().append('g').attr('class', 'root') enter.append('circle') .attr('r', 5) .attr('stroke-width', 0) .attr('stroke', 'transparent') .attr('fill', 'rgba(255, 255, 255, 0.7)') .attr('cx', d => currentXScale(d[0])) .attr('cy', d => currentYScale(d[1])) enter.append('circle') .attr('r', 3) .attr('stroke-width', 0) .attr('stroke', 'transparent') .attr('fill', 'rgba(0, 0, 0, 1)') enter.selectAll('circle') .attr('cx', d => currentXScale(d[0])) .attr('cy', d => currentYScale(d[1])) if (!debugRootFinder) { searchOutput = [{z: searchCircle.z, r: searchCircle.r, text: ''}] } let searchCircleSelection = gSearchCircles.selectAll('g.searchCircle') .data(searchOutput) .join(enter => { var appended = enter.append('g') .attr('class', 'searchCircle') appended.append('circle') .attr('stroke-width', 1) .attr('stroke', 'rgba(0, 0, 0, 0.8)') .attr('fill', 'none'); appended.append('text') .attr("font-style", "italic") .attr("font-size", "17px") .attr('dx', 10) .attr('alignment-baseline', 'middle') .attr('stroke', 'rgba(255, 255, 255, 0.7)') .attr('stroke-width', 3) .attr('paint-order', 'stroke') .attr('stroke-linecap', 'butt') .attr('stroke-linejoin', 'miter') .attr('stroke-miterlimit', 2); }) searchCircleSelection.exit().remove() gSearchCircles.selectAll('g.searchCircle') .select('circle') .attr('cx', d => currentXScale(d.z[0])) .attr('cy', d => currentYScale(d.z[1])) .attr('r', d => currentXScale(d.r) - currentXScale(0)) .attr('stroke', d => { return 'rgba(0, 0, 0, 1.0)' }); gSearchCircles.selectAll('g.searchCircle').select('text') .attr('x', d => currentXScale(d.z[0])) .attr('y', d => currentYScale(d.z[1])) .text(d => d.text) handles.selectAll('circle') .attr('cx', d => currentXScale(searchCircle.z[0] + (d === 'z0' ? 0 : searchCircle.r))) .attr('cy', d => currentYScale(searchCircle.z[1])) } function updateAxes () { axes.call(viewportAxes(viewport, currentXScale, currentYScale)) } const drawFunction = createDrawFunction(regl, fLambWave) function updatePlot () { updateRoots() } function updateFunction () { regl.poll() regl.clear({color: [1, 1, 1, 1]}) configureRegl.viewport(viewport, () => { configureRegl.scale(currentXScale, currentYScale, () => { configureRegl.map(() => { drawFunction({ω2, γ2}) }) }) }) } handles.selectAll('g.handles') .data(['z0', 'r0']) .enter().append('g') .attr('class', 'handle') .call(el => { el.append('circle') .attr('r', 6) .attr('stroke-width', 2) .attr('stroke', 'rgba(255, 255, 255, 0.9)') .attr('fill', '#35e') el.append('circle') .attr('cursor', 'move') .attr('r', 40) .attr('stroke-width', 0) .attr('stroke', 'transparent') .attr('fill', 'transparent') .call(d3.drag().on('drag', d => { if (d === 'r0') { searchCircle.r = currentXScale.invert(d3.event.x) - searchCircle.z[0]; } else { searchCircle.z[0] = currentXScale.invert(d3.event.x) searchCircle.z[1] = currentYScale.invert(d3.event.y) } computeRoots() updatePlot() })); }) svg.call(persistentZoom(currentXScale, currentYScale, xScale, yScale) .scaleExtent([1, 1]) .on('zoom.axes', updateAxes) .on('zoom.plot', updatePlot) .on('zoom.fn', updateFunction)) computeRoots() updateAxes(); updatePlot(); updateRoots() updateFunction() yield stack.container }3f:T522,md`## Putting it all together If we assemble the pieces we've learned so far—subdivision, contour integration to compute zeros, and Newton iteration to refine _pretty good_ guesses—we have something that starts to work. Our procedure looks like this: 1. Compute ${$`s_0`} for a given contour 2. If there are too many zeros or ${$`s_0`} is not close to an integer, subdivide the contour 3. For each subdivided contour, compute ${$`s_1, s_2, \ldots, s_Z`} 4. Use Newton's identities to construct a polynomial which shares the roots of the complex function 5. Pass the constructed polynomial to a standard polynomial root-finder. 6. Refine the computed roots with Newton's method. Naturally, this sequence of steps glosses over a number of little details like ensuring we don't accidentally compute the same zeros twice. Additionally, since we deal in large powers of ${$`z`}, it's better to integrate around the *unit* circle and transform variables to compute ${$`s_N`}. Finally, Ioakidimis et al. suggest avoiding the need for the function's derivative in their paper [A Modification of the Delves-Lyness Method for Locating the Zeros of Analytic Functions](https://www.sciencedirect.com/science/article/pii/00219991859012510). It sounds generally promising, but I haven't implemented that extension here. `40:T16a6,{ const { stack, regl, xScale, yScale, currentXScale, currentYScale, container, configureRegl, viewport, searchCircle } = delvesLynessPlot; const svg = d3.select(stack.layers.svg); svg.selectAll('g.plot, g.axes, defs').remove(); const plot = svg .append('g') .attr('class', 'plot') .attr('clip-path', 'url(#clip-rect)'); var searchOutput = []; var roots; function computeRoots() { roots = findComplexRoots(fLambWave.f, searchCircle.z, searchCircle.r, { integrationTolerance, maxIntegrationDepth, searchOutput, maxRootsPerCircle, maxSearchDepth }); if (refineRoots) { for (let i = 0; i < roots.length; i++) { complexNewtonRaphson(roots[i], fLambWave.f, roots[i]); } } } let axes = svg.append('g').attr('class', 'axes'); const gRoot = plot.append('g').attr('class', 'roots'); const gSearchCircles = plot.append('g').attr('class', 'searchCircles'); let handles = plot.append('g').attr('class', 'handles'); svg .append('defs') .append('clipPath') .attr('id', 'clip-rect') .append('rect') .attr('x', viewport.margin.l) .attr('y', viewport.margin.t) .attr('width', viewport.width - viewport.margin.l - viewport.margin.r) .attr('height', viewport.height - viewport.margin.t - viewport.margin.b); let searchCircles; function updateRoots() { gRoot.selectAll('g.root').remove(); const enter = gRoot .selectAll('g.root') .data(roots) .enter() .append('g') .attr('class', 'root'); enter .append('circle') .attr('r', 5) .attr('stroke-width', 0) .attr('stroke', 'transparent') .attr('fill', 'rgba(255, 255, 255, 0.7)') .attr('cx', d => currentXScale(d[0])) .attr('cy', d => currentYScale(d[1])); enter .append('circle') .attr('r', 3) .attr('stroke-width', 0) .attr('stroke', 'transparent') .attr('fill', 'rgba(0, 0, 0, 1)'); enter .selectAll('circle') .attr('cx', d => currentXScale(d[0])) .attr('cy', d => currentYScale(d[1])); if (!debugRootFinder) { searchOutput = [{ z: searchCircle.z, r: searchCircle.r, text: '' }]; } let searchCircleSelection = gSearchCircles .selectAll('g.searchCircle') .data(searchOutput) .join(enter => { var appended = enter.append('g').attr('class', 'searchCircle'); appended .append('circle') .attr('stroke-width', 1) .attr('stroke', 'rgba(0, 0, 0, 0.8)') .attr('fill', 'none'); appended .append('text') .attr("font-style", "italic") .attr("font-size", "17px") .attr('dx', 10) .attr('alignment-baseline', 'middle') .attr('stroke', 'rgba(255, 255, 255, 0.7)') .attr('stroke-width', 3) .attr('paint-order', 'stroke fill') .attr('stroke-linecap', 'butt') .attr('stroke-linejoin', 'miter') .attr('stroke-miterlimit', 2); }); searchCircleSelection.exit().remove(); gSearchCircles .selectAll('g.searchCircle') .select('circle') .attr('cx', d => currentXScale(d.z[0])) .attr('cy', d => currentYScale(d.z[1])) .attr('r', d => currentXScale(d.r) - currentXScale(0)) .attr('stroke', d => { var r = currentXScale(d.r) - currentXScale(0); return (r < 100 || r > 200) && debugRootFinder ? 'rgba(0, 0, 0, 0.2)' : 'rgba(0, 0, 0, 1.0)'; }); gSearchCircles .selectAll('g.searchCircle') .select('text') .attr('x', d => currentXScale(d.z[0])) .attr('y', d => currentYScale(d.z[1])) .text(d => { var r = currentXScale(d.r) - currentXScale(0); if (r < 100 || r > 200) return ''; return d.text; }); handles .selectAll('circle') .attr('cx', d => currentXScale(searchCircle.z[0] + (d === 'z0' ? 0 : searchCircle.r)) ) .attr('cy', d => currentYScale(searchCircle.z[1])); } function updateAxes() { axes.call(viewportAxes(viewport, currentXScale, currentYScale)); } const drawFunction = createDrawFunction(regl, fLambWave); function updatePlot() { updateRoots(); regl.poll(); regl.clear({ color: [1, 1, 1, 1] }); configureRegl.viewport(viewport, () => { configureRegl.scale(currentXScale, currentYScale, () => { configureRegl.map(() => { drawFunction({ ω2, γ2 }); }); }); }); } handles .selectAll('g.handles') .data(['z0', 'r0']) .enter() .append('g') .attr('class', 'handle') .call(el => { el.append('circle') .attr('r', 6) .attr('stroke-width', 2) .attr('stroke', 'rgba(255, 255, 255, 0.9)') .attr('fill', '#35e'); el.append('circle') .attr('cursor', 'move') .attr('r', 40) .attr('stroke-width', 0) .attr('stroke', 'transparent') .attr('fill', 'transparent') .call( d3.drag().on('drag', d => { if (d === 'r0') { searchCircle.r = currentXScale.invert(d3.event.x) - searchCircle.z[0]; } else { searchCircle.z[0] = currentXScale.invert(d3.event.x); searchCircle.z[1] = currentYScale.invert(d3.event.y); } computeRoots(); updatePlot(); }) ); }); svg.call( persistentZoom(currentXScale, currentYScale, xScale, yScale) .scaleExtent([0.5, 30]) .on('zoom.axes', updateAxes) .on('zoom.plot', updatePlot) ); computeRoots(); updateAxes(); updatePlot(); updateRoots(); yield stack.container; }41:T824,md`## Where do we go from here? So this is something. We found some roots without meandering about too haphazardly in the complex plane. Is the algorithm good? To be honest, it's not great. It's not the worst, but it's still inefficient compared to the state of the art, the inversion is, in general, ill-conditioned, and it struggles to distinguish what's actually present even for relatively straightforward cases (like poles or square roots). What could we improve? Well for one, we've integrated the function directly with a pretty standard [adaptive numerical integrator](https://observablehq.com/@rreusser/integration). In doing that, we've sampled points *non-uniformly* around the unit circle. What if we instead sampled the points *uniformly* around the unit circle? If you're familiar with the [Fast Fourier Transform](https://en.wikipedia.org/wiki/Fast_Fourier_transform) (FFT), talk of complex values uniformly spaced around the unit circle ("roots of unity") might perhaps ring a bell. Suddenly our analysis puts us FFT land, which is always a good sign. The bottom line is that it's a deep rabbit hole, and more sophisticated methods take advantage of additional structure to compute more valuable information with fewer function evaluations. [Jacob Rus](https://observablehq.com/@jrus) has directed me to some excellent papers including Austin et al.'s [Numerical algorithms based on analytic function values at roots of unity](http://people.maths.ox.ac.uk/trefethen/publication/PDF/2015_153.pdf) and Nakatsukasa et al.'s [The AAA algorithm for rational approximation](https://arxiv.org/pdf/1612.00337.pdf). They even include concise (though not permissibly licensed, I believe?) code to implement the methods. I really hope this has been interesting and has provided a small window into the intuition and mental picture behind some interesting math, but as far as writing numerical software goes, the rabbit hole is deep and I'll keep my day job for now. In the future, I'll use this technique to help explore the physical waves the pictured equation represents.`42:T8ec,md`## Addendum: The Lamb wave dispersion relation as a fun equation to work with We've been working with an abstract function, but it's not just a made up function. [Lamb waves](https://en.wikipedia.org/wiki/Lamb_waves) are planar elastodynamic waves that travel in a plate. For example, when you strike a metal plate and sound reverberates throughout—that's Lamb waves. Waves in unbounded media (like sound waves) are simple and may travel at any frequency. Waves in bounded media (like a metal plate) are significantly more complicated. The boundary conditions enforce constraints on which waves can travel so that only isolated spatial frequencies, or *modes*, may exist. A [dispersion relation](https://en.wikipedia.org/wiki/Dispersion_relation) relates spatial and temporal frequencies of traveling waves. The dispersion relation for sound in air is effectively just the speed of sound. For Lamb waves at a given temporal frequency, the dispersion relation describes the spatial frequencies of waves that can travel in the plate. The dispersion relations for symmetric (${$`+1`}) and antisymmetric (${$`-1`}) Lamb waves are given by ${$$`\displaystyle\frac{\tanh(\beta d / 2)}{\tanh(\alpha d / 2)} = \left( \frac{4 \alpha \beta k^2}{(k^2 + \beta^2)^2} \right)^{\pm 1}`} with ${$$`\displaystyle \alpha^2 = k^2 - \frac{\omega^2}{c_l^2},`} ${$$`\beta^2 = k^2 - \frac{\omega^2}{c_t^2}`} where ${tex`d`} is the thickness of the plate, ${tex`k`} the wavenumber in the plane of the plate, ${tex`\omega`} the temporal frequency, ${tex`c_l`} the longitudinal wave speed, and ${tex`c_t`} the shear (transverse) wave speed. As noted above, we of course don't actually implement this function directly but cross-multiply to remove any possible poles. The branch cuts resulting from the implicit square roots in ${$`\alpha`} and ${$`\beta`} are somewhat more problematic. Squaring the entire equation would resolve the branch cuts at the cost of repeating all other roots, so we grin and bear the branch cuts. Wavenumbers ${tex`k`} which satisfy this equation represent the different superpositions of bulk waves that satisfy the boundary conditions of the plate. That is, in order to determine which waves are valid, we must find the zeros of this equation in the complex plane.`43:T908,fLambWave = ({ glsl: ` // This function puts the branch cut on the other side of the square // root, away from most of the function we're looking at. It's // mathematically equivalent to csqrt though. complex csqrtAlt (complex a) { return cmul(vec2(0, 1), csqrt(-a)); } // We write this function generically and replace 'complex' with either // vec2 or vec4, depending on whether we're using automatic differentiation complex f (complex k2) { k2 = csqr(k2); complex alpha2 = csub(k2, w2); complex beta2 = csub(k2, cmul(w2, gamma2)); complex alpha = csqrtAlt(alpha2); complex beta = csqrtAlt(beta2); complex eb = cexp(beta); complex enb = cexp(-beta); complex ea = cexp(alpha); complex ena = cexp(-alpha); // Division by two cancels out, so we omit it complex sinhb2 = (eb - enb); complex coshb2 = (eb + enb); complex sinha2 = (ea - ena); complex cosha2 = (ea + ena); complex k2b22 = csqr(cadd(k2, beta2)); complex ab4k2 = 4.0 * cmul(cmul(alpha, beta), k2); complex lhs = cmul(cmul(sinhb2, cosha2), ${parity === 'Symmetric' ? 'k2b22' : '-ab4k2'}); complex rhs = cmul(cmul(sinha2, coshb2), ${parity === 'Symmetric' ? 'ab4k2' : '-k2b22'}); return rhs - lhs; } `, createUniforms: function (regl) { return { w2: regl.prop('ω2'), gamma2: regl.prop('γ2'), } }, f: (function () { var tmp = ccreateD() var k2 = ccreateD() var alpha2 = ccreateD(), beta2 = ccreateD() var alpha = ccreateD(), beta = ccreateD() var sinha = ccreateD(), cosha = ccreateD() var sinhb = ccreateD(), coshb = ccreateD() var ab4k2 = ccreateD(), k2b22 = ccreateD() var lhs = ccreateD(), rhs = ccreateD() var symm = parity === 'Symmetric' return function fLamb (out, a, b) { csetValuesD(k2, a, b, 1, 0) csqrD(k2, k2) csubDS(alpha2, k2, ω2) csubDS(beta2, k2, cmulSS(tmp, ω2, γ2)) csqrtD(alpha, alpha2) csqrtD(beta, beta2) csinhcoshD(sinha, cosha, alpha) csinhcoshD(sinhb, coshb, beta) csqrD(k2b22, caddDD(tmp, k2, beta2)) cscaleD(ab4k2, cmulDD(tmp, cmulDD(tmp, alpha, beta), k2), 4.0) cmulDD(lhs, cmulDD(tmp, sinhb, cosha), symm ? k2b22 : ab4k2) cmulDD(rhs, cmulDD(tmp, sinha, coshb), symm ? ab4k2 : k2b22) return csubDD(out, lhs, rhs) } })() })44:T507,function createPlot (config) { var xrange = config.xrange || [-10, 10] var yrange = config.yrange || [-10, 10] var margins = config.margins || {t: 10, r: 5, b: 20, l: 30} const stack = createLayerStack(width, Math.max(500, Math.min(900, width * 0.7)), Math.min(devicePixelRatio), { regl: reglCanvasWithOptions({ attributes: {depthStencil: false, antialias: false}, optionalExtensions: ['OES_standard_derivatives'], }), svg: DOM.svg }) const viewport = createViewport(stack, margins) const yScale = d3.scaleLinear() .domain(xrange) .range([viewport.height - viewport.margin.b, viewport.margin.t]) const xScale = constrainLinearScaleAspectRatio( d3.scaleLinear() .domain(yrange) .range([viewport.margin.l, viewport.width - viewport.margin.r]), yScale, 1) const configureRegl = ({ viewport: createReglViewportConfiguration(stack.layers.regl), scale: createReglLinearScaleConfiguration(stack.layers.regl), map: createReglMap(stack.layers.regl), }) return Object.assign({ stack, container: stack.container, regl: stack.layers.regl, viewport, xScale, yScale, currentXScale: xScale.copy(), currentYScale: yScale.copy(), configureRegl, }, config.extras) }45:T7b2,function createDrawFunction(regl, func) { return regl({ frag: ` ${ derivatives === "OES_standard_derivatives" ? "#extension GL_OES_standard_derivatives : enable" : "" } precision highp float; uniform vec2 gamma2, w2; uniform sampler2D colormap; varying vec2 xy; uniform vec2 inverseViewportResolution; uniform mat4 inverseView; uniform float n; ${createPolarDomainColoringShader({ magnitudeOctaves: 4, phaseOctaves: 4, argumentColoring: 1.0, bias1: 0.5 })} ${glslComplexAutoDifferentiation} ${glslRainbow} // All of the functions are generalized to work with either vec2 (re, im) // or with vec4 (the function and its derivative). We generalize the type // and nail two birds with one stone. ${func.glsl.replace( /\bcomplex\b/g, derivatives === "OES_standard_derivatives" ? "vec2" : "vec4" )} void main () { ${ derivatives === "OES_standard_derivatives" ? ` vec2 fz = f(xy); vec4 fdf = vec4(fz, fwidth(fz) * 0.5); ` : ` float unitsPerPixel = inverseView[0][0] * inverseViewportResolution.x * 2.0; vec4 fdf = f(vec4(vec2(xy), vec2(unitsPerPixel, 0))); ` } gl_FragColor = mix(vec4(vec3(1.0), 1.0), domainColoring( fdf, vec2(6), // steps vec2(0.2), // scale vec2(0.2, 0.3), // grid vec2(0.5, 0.1), // shading 1.0, 1.0, // line width, line feather vec3(0.0), // grid color 1.0, // argument coloring 1.0, // contrast power colormap, vec2(0) ), 0.9); }`, uniforms: Object.assign( func.createUniforms ? func.createUniforms(regl) : {}, { colormap: createRainbowTexture(regl) } ) }); }46:T69f,function createFractalDrawFunction(regl, func) { return regl({ frag: ` precision highp float; varying vec2 xy; uniform vec2 inverseViewportResolution; uniform mat4 inverseView; uniform sampler2D colormap; ${createPolarDomainColoringShader({ magnitudeOctaves: 4, phaseOctaves: 4, argumentColoring: 0.001, bias1: 0.5 })} ${glslComplexAutoDifferentiation} ${glslRainbow} ${func.glsl.replace(/\bcomplex\b/g, "vec4")} void main () { vec2 z = xy; float lightness = 1.0; for (int i = 0; i < 40; i++) { vec4 fdf = f(vec4(z, 1.0, 0.0)); z = z - cdiv(fdf.xy, fdf.zw); lightness *= mix(1.0, 0.988, smoothstep(-2.5, 1.0, log(dot(fdf.xy, fdf.xy)))); } lightness *= lightness; lightness *= lightness; float unitsPerPixel = inverseView[0][0] * inverseViewportResolution.x * 2.0; vec4 fdf = f(vec4(vec2(xy), vec2(unitsPerPixel, 0))); vec3 dc = domainColoring( fdf, vec2(6), // steps vec2(0.2), // scale vec2(0.2, 0.2), // grid vec2(0.15, 0.0), // shading 1.0, 1.0, // line width, line feather vec3(0.0), // grid color 0.0, // argument coloring 1.0, // contrast power colormap, vec2(0) ).xyz; gl_FragColor = vec4(rainbow(atan(z.y, z.x) * 0.5 / 3.14159 + 0.25) * lightness * dc, 1.0); }`, uniforms: Object.assign( func.createUniforms ? func.createUniforms(regl) : {}, { colormap: createRainbowTexture(regl) } ) }); }47:T4a1,md`## License With the exception of code imported from other sources as stated above, the code in this notebook is MIT Licensed. Copyright 2019 Ricky Reusser Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. `48:[] 30:{"id":"6328fc7ba0803522","slug":"finding-roots-in-the-complex-plane","title":"Finding Roots in the Complex Plane","thumbnail":"277c16eb64b16b63b2b8754c287ca27729a6d6d1f33ae1260de47b178436d117","access_level":"public","create_time":"$D2019-07-27T02:34:52.298Z","update_time":"$D2023-09-05T04:22:23.317Z","view_time":null,"first_public_time":"$D2019-08-19T06:15:02.021Z","has_unpublished_edits":null,"fork":false,"stars":59,"starred":false,"owner":{"type":"team","id":"f1b8352b963a77b4","login":"rreusser","name":"Ricky Reusser","avatar_url":"https://avatars-next.observableusercontent.com/49c18a2ac0092176d4ea8348a1dd538e0b460b086b89d602b8d21ffe0821d92e","tier":"starter_2024"},"version":6160,"publish_version":null,"role":null,"authored":false,"authors":[{"type":"individual","id":"908c4be2279a714e","login":"rreusser","name":"Ricky Reusser","avatar_url":"https://avatars-next.observableusercontent.com/fa41989f47bb986723af5b9c22176cb503748426f56e565000ea93acd9eab86c","tier":null}],"members":[{"type":"individual","id":"908c4be2279a714e","login":"rreusser","name":"Ricky Reusser","avatar_url":"https://avatars-next.observableusercontent.com/fa41989f47bb986723af5b9c22176cb503748426f56e565000ea93acd9eab86c","tier":null,"role":"editor"}],"comments":[],"custom_thumbnail":true,"trashed":false,"lapsed":null,"subscribed":null,"forks":10,"collection":{"id":"85cbbc7776cf19d5","slug":"math","title":"Math","owner":{"id":"f1b8352b963a77b4","login":"rreusser","name":"Ricky Reusser"}},"type":"notebook","body":{"titling":"auto","title":"Finding Roots in the Complex Plane","cells":[{"id":3271,"value":"$32","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"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

🌈 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","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"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","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"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","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"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","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":4831,"value":"$36","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"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","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"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","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"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","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":0,"value":"$37","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"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","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":3477,"value":"$38","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"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","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":3657,"value":"md`## Problems and annoyances`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"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","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":3670,"value":"$39","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"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","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"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","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"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","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":4152,"value":"$3a","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"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","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":3865,"value":"$3b","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":4062,"value":"$3c","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"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","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":1140,"value":"$3d","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":3970,"value":"$3e","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"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","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":4509,"value":"$3f","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":3608,"value":"$40","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"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","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":2010,"value":"viewof debugRootFinder = checkbox({options: [{label: 'Visualize search hierarchy', value: 'debug'}], value: 'debug'})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":2946,"value":"viewof refineRoots = checkbox({options: [{label: 'Refine roots with Newton-Raphson', value: 'refine'}], value: 'refine'})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"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","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":2741,"value":"viewof maxSearchDepth = slider({min: 1, max: 16, step: 1, value: 7, description: 'maximum search depth'})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"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","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"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","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":4774,"value":"$41","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"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","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":3157,"value":"$42","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":91,"value":"viewof f = slider({min: 100, max: 2e6, step: 100, value: 500e3, description: 'frequency, f, Hz'})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"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","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"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","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"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","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":5286,"value":"md`Plate thickness, ${$`d`}, in meters`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":13,"value":"d = 10e-3 // meters","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":5291,"value":"md`Density, ${$`\\rho`} in ${$`kg/m^3`}`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":9,"value":"ρ = 2700 // kg / m^3","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":5323,"value":"md`Magnitude of Young's modulus, ${$`E`}, in ${$`Pa`}`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":963,"value":"Emag = 68.9e9 // Pascals","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"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","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":1346,"value":"E = cscaleS([], [Math.cos(θ), Math.sin(θ)], Emag)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"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","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":1360,"value":"cl = cscaleS([], csqrtS([], E), Math.sqrt((1 - ν) / ρ / (1 + ν) / (1 - 2 * ν)))","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":5334,"value":"md`Shear (tangential) wave speed, ${$`c_t`}, in ${$`m/s`}.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":1367,"value":"ct = cscaleS([], csqrtS([], E), 1 / Math.sqrt(2 * ρ * (1 + ν)))","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":1387,"value":"ω2 = csqrS([], cscaleS([], cinvS([], cl), 2.0 * Math.PI * f * d)) ","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":1394,"value":"γ2 = csqrS([], cdivSS([], cl, ct))","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":101,"value":"$43","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":4406,"value":"md`## Plotting helpers`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":3920,"value":"$44","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":72,"value":"$45","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":6113,"value":"import { rainbow } from '481d4902e37d39c8'","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":6115,"value":"createRainbowTexture = regl => regl.texture(rainbow)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":4833,"value":"$46","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"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","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"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","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":76,"value":"import { glslRainbow } from '@rreusser/adaptive-domain-coloring'","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":6109,"value":"import { createPolarDomainColoringShader } from '39fe342decb60c30'","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":4396,"value":"md`## Math helpers`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"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","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"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","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":3109,"value":"import {adaptiveSimpson, vectorAdaptiveSimpson} from '@rreusser/integration'","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":1187,"value":"import { findRoots } from '@rreusser/polynomial-roots'","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":78,"value":"import {glslComplexAutoDifferentiation} from '@rreusser/glsl-complex-auto-differentiation'","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":50,"value":"d3 = require('d3@5')","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":89,"value":"import {slider, radio, checkbox} from '@jashkenas/inputs'","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":4451,"value":"$$ = tex","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":4448,"value":"$$$ = $.block","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":3597,"value":"html``","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":3822,"value":"$47","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"}],"files":"$W48","stdlib":"1"}} 31:"$30:body"