Inversion in a circle

Inversion in the unit circle is a way of turning the circle inside-out. Everything that was inside the circle goes outside the circle, and everything that was outside the circle comes in.

Not only is the disk turned inside-out, the same thing happens along each ray going out from the origin. Points on that ray that are inside the circle go outside and vice versa. In polar coordinates, the point (r, θ) goes to (1/r, θ).

Complex numbers

In terms of complex numbers, inversion in the unit circle amounts to taking the reciprocal and the conjugate (in either order, because these operations commute). This is the same as dividing a complex number by the square of its magnitude. Proof:

z \bar{z} = |z|^2 \implies \frac{1}{\bar{z}} = \frac{z}{|z|^2}

There are two ways to deal with the case z = 0. One is to exclude it, and the other is to say that it maps to the point at infinity. This can be made rigorous by working on the Riemann sphere rather than the complex plane. More on that here.

Inverting a hyperbola

The other day Alex Kantorovich pointed out on Twitter that “the perfect figure 8 (or infinity) is simply the inversion through a circle of a hyperbola.” We’ll demonstrate this with Mathematica.

What happens to a point on the hyperbola

x^2 - y^2 = 1

when you invert it through a circle? If we think of x and y as the real and imaginary parts of a complex number, the discussion above shows that the point (x, y) goes to the same point divided by its length squared.

(x, y) \mapsto \left( \frac{x}{x^2 + y^2}, \frac{y}{x^2 + y^2} \right)

Let (u, v) be the image of the point (x, y).

\begin{align*} u &= \frac{x}{x^2 + y^2} \\ v &= \frac{y}{x^2 + y^2} \end{align*}

Then we have

u^2 + y^2 = \frac{x^2 + v^2}{(x^2 + y^2)^2} = \frac{1}{x^2 + y^2}

and

u^2 - v^2 = \frac{x^2 - y^2}{(x^2 + y^2)^2} = \frac{1}{(x^2 + y^2)^2}

because x² − y² = 1. Notice that the latter is the square of the former, i.e.

u^2 - v^2 = (u^2 + v^2)^2

Now we have everything to make our plot. The code

ContourPlot[{
    x^2 + y^2 == 1, 
    x^2 - y^2 == 1, 
    x^2 - y^2 == (x^2 + y^2)^2}, 
    {x, -3, 3}, {y, -3, 3}]

produces the following plot.

The blue circle is the first contour, the orange hyperbola is the second contour, and the green figure eight is the third contour.

Related posts

Visual integration

The plot below is of a meromorphic function f(z). That is, the function f(z) is analytic except possibly at poles, and the colors represent the phase angles, the values of θ if you write the function values in polar form.

What is the value of the integral

\frac{1}{2\pi i} \int_C \frac{f'(z)}{f(z)}\, dz

where C is the perimeter of the square?

This hardly seems like enough information to go on. I haven’t even said what the function is. And yet all the information you need is there.

There are three important points inside the plot, the three points where the color swirls come to a point. The colors rotate around the point near the top right corner twice, and they rotate around the other two points once each, and in the opposite order.

This means that our function either has two zeros and a double pole, or two poles and a double zero. In order to know which is which we need to know the color convention. The plot was made with Mathematica, and colors go through ROYGBIV order counterclockwise as the angle increases. So we have the former case: two single zeros and one double pole. However, it turns out the integral would be the same in either case.

By the argument principle, the integral above counts the number of zeros minus the number of poles inside the region, counted with multiplicity. Each zero contributes +1 and the double pole contributes −2. The integral is zero.

Now suppose we cut the region along a diagonal from the top left to the bottom right. Let C be the perimeter of the right triangle below the diagonal and integrate around C counterclockwise. What’s the value of the integral?

It would be 2, because there are two zeros inside the triangle.

***

To learn more about interpreting contour plots, see this post. For more on the argument principle, read this post then this post.

In my previous posts on the argument principle, I only mentioned zeros. Everything holds for poles as well, if you flip the sign.

Sum the zeros of an analytic function without finding them first

A couple days ago I wrote about how Vieta’s formulas let you sum the zeros of a polynomial without having to first compute the zeros. This is especially handy for high-order polynomials since there is no explicit formula for the zeros.

Most functions that arise in applications are not polynomials. How could you find the sum of the zeros of an analytic function in some region without having to locate each of the zeros first? What if you don’t even know how many zeros are in the region?

Here’s a plot of the complex function we’re going to use as an example.

Complex 3D plot of Bessel function Y_1

As with yesterday’s post, the key is to use the argument principle, but this time a more general form.

Generalized argument principle

Suppose we have a simple closed curve C and a function f(z) that is analytic in and on C. In particular we’re assuming f has no poles in or on C. Assume f has n zeros labeled z1 through zn. We’re going to integrate around C to add up the zeros of f inside.

Let g(z) be a function that we want to apply to the zeros of f. We assume it is also analytic in and on C. Then

\frac{1}{2\pi i}\int_C g(z)\, \frac{f'(z)}{f(z)}\, dx = \sum_{i=1}^n g(z_i)

In words, we can sum the values of g applied to each of the zeros of f by computing the integral on the left.

In the previous post, we had g(z) = 1. That is, we were simply counting the zeros. To sum the zeros, we set g(z) = z.

Example

The function plotted above is the Bessel function Y1. Here’s a flattened plot, just showing the phase of the values.

Contour plot of Bessel function Y_1

The points along the real line where the colors swirl around a point are the zeros of the function. We can see that there are zeros near 2, 5, 9, and 11.

There’s something different going on at 0. The colors swirl in the opposite direction because there’s a pole at 0. And there’s an abrupt change in color on the real line to the left of 0 because there’s a branch cut there.

We want to find the sum of the first four zeros, so we’ll create a contour that includes these zeros and excludes the singularity at 0. We’ll use a rectangle with lower left corner at 1 − i and an upper right corner at 13 + i.

We’ll modify our Mathematica code from yesterday to use a different example function, Y1, and add a factor of z to the integrand to sum zeros rather than count zeros.

    f[z_] := BesselY[1, z]
    leg[z1_, z2_] := NIntegrate[z f'[z]/f[z], {z, z1, z2}]
    (leg[1-I, 13-I] + leg[13-I, 13+I] + leg[13+I, 1+I] + leg[1+I, 1-I])/(2 Pi I)

This returns 27.972 − 5.65432*10^−16 I. We know that the zeros are all real, so the imaginary part is integration error. We can ask Mathematica for more digits, and it will report the real part as 27.97198306597801.

We can compute the sum of the first four zeros directly

    N[Sum[BesselYZero[1, n], {n, 1, 4}]]

and the result agrees with the result above to 13 significant figures.

When we were counting zeros, we knew the result had to be an integer, so the rounded integral was exact. Here our result is only as accurate as our numerical integration, though we did use prior information to know that we could discard the imaginary part of the integration result. This would not have made any difference since the imaginary part is orders of magnitude smaller than the integration error in the real part.

Related posts

Bounding zeros of an analytic function

The previous post looked at the problem of finding the zeros of a cubic polynomial. Assuming we’re going to use a numerical method to calculate the zero, the hard part is knowing where to tell the numerical method to look. That post showed how to use a change of variables to guarantee that the polynomial has one, and only one, root in the unit interval.

Now suppose we’re looking at functions of a complex variable. On the real line we can say that if a function is negative here and positive there, it must be zero somewhere in between by the intermediate value theorem. But that’s not the case for a complex-valued function. Such a function can avoid going through zero because it can move in two dimensions, not just one.

For example, the function ez is never zero. Now eπi = −1 and e0 = 1, but there’s no point on a line from πi to 0 where the exponential function is zero.

So how can we tell whether an analytic function has a zero in some region of the complex plane? More generally, can we tell how many zeros it has in some region?

There is a theorem in complex analysis called the argument principle. It says that if an analytic function f is not zero along a closed curve C, then the number of zeros of f inside C equals

\frac{1}{2\pi i} \int_C \frac{f'(z)}{f(z)}\, dz

We can evaluate this integral numerically, and it will tell us the number of zeros inside. The exact value must be an integer, so the integral doesn’t have to be computed with much accuracy. If we know the error is less than a half, then rounding the result to the nearest integer gives the exact answer.

Riemann zeta example

Let’s use the argument principle to show that the Riemann zeta function ζ(z) has 3 zeros in the critical strip with imaginary part between 10 and 30. The critical strip is the part of the complex plane with real part between 0 and 1.

We can make our contour a rectangle with lower left corner at 10i and upper right corner at 1 + 30i.

We’ll use Mathematica and start by defining a function leg that numerically integrates along one leg of a rectangle.

    leg[z1_, z2_] := NIntegrate[Zeta'[z]/Zeta[z], {z, z1, z2}]

Now let’s use this function to integrate over the rectangle described above.

    (leg[0+10I, 1+10I] + leg[1+10I, 1+30I] + 
     leg[1+30I, 0+30I] + leg[0+30I, 0+10I]) / (2 Pi I)

This returns

    3. +3.77069*10^-12 I

and the nearest integer is 3, confirming that the zeta function has 3 zeros inside our box, assuming that our numerical integration error is less than 0.5.

Next, let’s show that zeta has a zero in the critical strip with imaginary part between 14 and 15.

The integration

    (leg[0+14I, 1+14I] + leg[1+14I, 1+15I] + 
     leg[1+15I, 0+15I] + leg[0+15I, 0+14I]) / (2 Pi I)

returns

    1. + 1.46221*10^-12 I

showing that there’s one zero inside this smaller box.

In fact the first three zeros in the critical strip are 1/2 + 14.1347i, 1/2 + 21.0220i, and 1/2 + 25.0109i. The next one is at 1/2 + 30.4249i, just outside the first box we used above.

Here’s a plot.

I plotted ζ(1/2 + iz) to turn the plot sideways. Here’s the code that produced the plot.

    ComplexPlot[Zeta[I z + 1/2], {z, 10 - I, 30 + I}, 
        ColorFunction -> "LocalMaxAbs", AspectRatio -> 1/3]

More on plots of a complex function and what the colors mean here.

Related posts

Reading plots of a complex function

This post is about understanding the following plot. If you’d like to read more along these lines, see [1].

Hankel function phase plot

The plot was made with the following Mathematica command.

    ComplexPlot[HankelH1[3, z], {z, -8 - 4 I, 4 + 4 I}, 
        ColorFunction -> "CyclicArg", AspectRatio -> 1/2]

The plot uses color to represent the phase of the function values. If the output of the function is written in polar form as (r, θ), we’re seeing θ.

Here’s a plot of the identity function f(z) = z to show how the color rotation works.

The color rotate from red to orange to green etc. (ROYGBIV) clockwise around a pole and counterclockwise around a zero.

You can see from the plot at the top that our function has a pole at 0. In fact, it’s a pole of order three because you go through the ROYGBIV cycle three times clockwise as you circle the origin.

You can also see four zeros, located near −6.4 − 0.4i, −2, 2 − i, −0.4 − 2i, and −1.3 − 1.7i. In each case we cycle through ROYGBIV one time counterclockwise, so these are simple zeros. If we had a double zero, we’d cycle through the colors twice. If we had a triple zero, we’d cycle through the colors three times, just as we did at the pole, except we’d go through the colors counterclockwise.

The colors in our plot vary continuously, except on the left side. There’s a discontinuity in our colors above and below the real axis. That’s because the function we’re plotting has a branch cut from −∞ to 0. The discontinuity isn’t noticeable near the origin, but it becomes more noticeable as you move away from the origin to the left.

Here’s a 3D plot to let us see the branch cut more clearly.

Here color represents phase as before, but now height represents absolute value, the r in our polar representation. The flat plateau on top is an artifact. The function becomes infinite at 0, so it had to be cut off.

You can see a seam in the graph. There’s a jump discontinuity along the negative real axis, with the function taking different values as you approach the real axis from the second quadrant or from the third quadrant.

This came up in a previous post, Architecture and math, but I didn’t say anything about it. Here’s a plot from that post.

Why is there a little radial line on top where the graph was chopped off? And why is there a white streak below it? That’s evidence of a branch cut, though the discontinuity is small in this graph. You’ll notice color swirls around the base of the mesa. The colors swirl counterclockwise because they’re all at zeros. But you’ll notice the colors swirl clockwise as you go around the walls of the mesa. In fact they swirl around 5 times because there’s a pole of order 5 at the origin.

The book I mentioned in that post had something like this on the cover, plotting a generalization of the Hankel function. That’s why I chose the Hankel function as my example in this post.

Related posts

[1] Elias Wegert. Visual Complex Functions: An Introduction with Phase Portraits

The ring of entire functions

Rings made a bad first impression on me. I couldn’t remember the definitions of all the different kinds of rings, much less have an intuition for what was important about each one. As I recall, all the examples of rings in our course were variations on the integers, often artificial variations.

Entire functions

I’m more interested in analysis than algebra, so my curiosity was piqued when I ran across an appendix on entire functions in the back of an algebra book [1]. This appendix opens with the statement

The ring E of entire functions is a most peculiar ring.

It’s interesting that something so natural from the perspective of analysis is considered peculiar from the perspective of algebra.

An entire function is a function of a complex variable that is analytic in the entire complex plane. Entire functions are functions that have a Taylor series with infinite radius of convergence.

Rings

A ring is a set of things with addition and multiplication operations, and these operations interact as you’d expect via distributive rules. You can add, subtract, and multiply, but not divide: addition is invertible but multiplication is not in general. Clearly the sum or product of entire functions is an entire function. But the reciprocal of an entire function is not an entire function because the reciprocal has poles where the original function has zeros.

So why is the ring of analytic functions peculiar to an algebraist? Osborne speaks of “the Jekyll-Hyde nature of E,” meaning that E is easy to work with in some ways but not in others. If Santa were an algebraist, he would say E is both naughty and nice.

Nice properties

On the nice side, E is an integral domain. That is, if f and g are entire functions and fg = 0, then either f = 0 or g = 0.

If we were looking at functions that were merely continuous, it would be possible for f to be zero in some places and g to be zero in the rest, so that the product fg is zero everywhere.

But analytic functions are much less flexible than continuous functions. If an analytic function is zero on a set of points with a limit point, it’s zero everywhere. If every point in the complex plane is either a zero of f or a zero of g, one of these sets of zeros must have a limit point.

Another nice property of E is that it is a Bezóut domain. This means that if f and g are entire functions with no shared zeros, there exist entire functions λ and μ such that

λf + μg = 1.

This is definition is analogous to (and motivated by) Bezóut’s theorem in number theory which says that if a and b are relatively prime integers, then there are integers m and n such that

ma + nb = 1.

Naughty properties

The naughty properties of E take longer to describe and involve dimension. “Nice” rings have small Krull dimension. For example Artinian rings have Krull dimension 0, and the ring of polynomials in n complex variables has dimension n. But the Krull dimension of the ring of entire functions is infinite. In fact it’s “very infinite” in the sense that it is at least

2^{\aleph_1}

and so the Krull dimension of E is larger than the cardinality of the complex numbers.

Wittgenstein’s ruler

Nassim Taleb described Wittgenstein’s ruler this way:

Unless you have confidence in the ruler’s reliability, if you use a ruler to measure a table you may also be using the table to measure the ruler. The less you trust the ruler’s reliability, the more information you are getting about the ruler and the less about the table.

An algebraist would say that entire functions are weird, but an analyst could say that on the contrary, ring theory, or at least Krull dimension, is weird.

Related posts

[1] Basic Homological Algebra by M. Scott Osborne.

Numbering the branches of the Lambert W function

The previous post used the Lambert W function to solve an equation that came up in counting partitions of an integer. The first solution found didn’t make sense in context, but another solution, one on a different branch, did. The default branch k = 0 wasn’t what we were after, but the branch k = -1 was.

Branches 0 and -1

The branch corresponding to k = 0 is the principal branch. In mathematical notation, the k is usually a subscript, i.e W0(z) is the principal branch of the Lambert W function. For real z ≥ -1/e it returns a real value. It’s often what you want, and that’s why it’s the default in software implementations such as SciPy and Mathematica. More on that below.

The only other branch that returns real values for real inputs is W-1 which returns real values for -1/ez < 0.

If you’re working strictly with real numbers, you only need to be concerned with branches 0 and -1. If branch 0 doesn’t give you what you expect, try branch -1, if your argument is negative.

SciPy and Mathematica

SciPy implements Wk with lambertw(z). The parameter k defaults to 0.

The Mathematica function ProductLog[z] implements W0(z), and ProductLog[k, z] implements Wk.

Note that Mathematica and Python list their arguments in opposite orders.

Branch numbering

The branches are numbered in order of their imaginary parts. That is, The imaginary part of Wk (z) is an increasing function of k. For example, here is a list of the numerical values of the imaginary parts of Wk (1) for k running from -3 to 3.

    Table[N[Im[ProductLog[k, 1]]], {k, -3, 3}]
    {-17.1135, -10.7763, -4.37519, 0., 4.37519, 10.7763, 17.1135}

Note the zero in the middle because W0(1) is real.

Recovering k

Given z and a value Wk (z) with k unknown, you can determine k with

k = \frac{W_k(z) + \log W_k(z) - \log z}{2\pi i}

with the exception that if the expression above is 0 and -1/ez < 0 then k = -1. See [1].

For example, let z = 1 + 2i and w = W3 (z).

    z = 1 + 2 I
    w = ProductLog[3, z]

Now pretend you don’t know how w was computed. When you compute

    N[(w + Log[w] - Log[z])/(2 Pi I)]

the result is 3.

Partition example

The discussion of the W function here grew out of a problem with partitions. Specifically, we wanted to solve for n such that
f(n) = \frac{1}{4n\sqrt{3}} \exp\left( \pi \sqrt{\frac{2n}{3}} \right)

is approximately a trillion. We found in the previous post that solutions of the equation

a w^b \exp(c w^d) = x

are given by

w = \left(\frac{b}{cd} \,\,W\left(\frac{cd}{b} \left(\frac{x}{a} \right )^{d/b} \right )\right )^{1/d}

Our partition problem corresponds to a = 1/√48, b = -1, c = π √(2/3), and d = 1/2. We found that the solution we were after came from the k = -1 branch.

    N[( (b/(c d)) ProductLog[-1, (x/a)^(d/b) c d /b])^(1/d)]
    183.852

Even though our final value was real, the intermediate values were not. The invocation of W-1 in the middle returns a complex value:

    N[ProductLog[-1, ((x/a)^(d/b) c d /b)^(1/d)]]
    -32.5568 - 3.24081 I

We weren’t careful about specifying ranges and branch cuts, but just sorta bulldozed our way to a solution. But it worked: it’s easy to verify that 183.852 is the solution we were after.

***

[1] Unwinding the branches of the Lambert W function by Jeffrey, Hare, and Corless. The Mathematical Scientist 21, 1&ndash;7 (1996)

Marden’s amazing theorem

The previous post was a warmup for this post. It gave an example of the theorem that if p is a polynomial, the roots of its derivative p′ lie inside the convex hull of the roots of p. If p is a cubic polynomial, we can say much more.

Suppose p(z) is a polynomial with three distinct roots, not all on a line. Then the roots of p are the vertices of a triangle T in the complex plane. We know that the roots of p′ lie inside T, but Marden’s theorem says where they lie inside T.

Let E be the unique ellipse inscribed inside T that touches T at the midpoint of each side. Then Marden’s theorem says that the roots of p′ lie on the foci of E.

The ellipse E is sometimes called the Steiner inellipse.

For an example, let

p(z) = z (z – 3) (z – 2 – 4i)

The roots of the derivative of p are 1.4495 + 0.3100i and 1.8838 + 2.3566i.

Here’s a plot of the roots of p (blue dots), the roots of p′ (orange ×), and the Steiner ellipse. I used a parametric equation for the Steiner ellipse from here.

Related posts

Convex hull of zeros

There’s a well-known theorem in complex analysis that says that if p is a polynomial, then the zeros of its derivative p′ lie inside the convex hull of the zeros of p. The convex hull of a set of points is the smallest convex set containing those points.

This post gives a brief illustration of the theorem. I created a polynomial with roots at 0, i, 2 + i, 3-i, and 1.5+0.5i. The convex hull of these points is the quadrilateral with corners at the first four roots; the fifth root is inside the convex hull of the other roots.

The roots are plotted with blue dots. The roots of the derivative are plotted with orange ×’s.

In the special case of cubic polynomials, we can say a lot more about where the roots of the derivative lie. That is the topic of the next post.

More complex analysis posts