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) ( − 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

Complex floor and a surprising pattern

The floor of a real number x, written ⌊x⌋, is the largest integer less than or equal to x. So, for example, ⌊π⌋ = 3 and ⌊−π⌋ = -4.

The previous post applied the floor function to a complex number z. What does that mean? You can’t just say it’s the largest integer [1] less than z because the complex numbers aren’t ordered; there’s no natural way to say whether one complex number is less than or greater than another.

The post was using Mathematica code, and Mathematica defines the floor of a complex number by applying the function to the real and imaginary parts separately. That is, if x and y are real numbers, then

x + iy⌋ = ⌊x⌋ + iy⌋.

That lets you take the floor of a complex number, but is it useful? I wondered how many identities that hold for floors of real numbers extend to complex numbers when the floor is defined as above, so I looked at the chapter in Concrete Mathematics that has a bunch of identities involving floors and ceilings to see which ones extend to complex numbers.

Any identity where the real an imaginary parts don’t interact should extend trivially. For example, for real x and integer n,

x + n⌋ = ⌊x⌋ + n.

This extends to the case where x is complex and n is a (complex) integer [1]. Addition works on real and imaginary components separately, just like Mathematica’s definition of floor.

But here’s one that’s much more interesting. For positive x,

⌊√⌊x⌋ ⌋ = ⌊√x⌋.

Taking square roots deeply intertwines the real and imaginary parts of a number [3]. Does the identity above hold for complex inputs?

I wrote some Python code to test this, and to my great surprise, the identity often holds. In my first experiment, it held something like 84% of the time for random inputs. I expected it would either rarely hold, or hold say half the time (e.g. in a half-plane).

My first thought was to plot 1,000 random points, green dots where the identity holds and red stars where it does not. This produced the following image.

Red stars scattered randomly among green dots

About 16% of the points are red and the rest green. In this plot there’s no apparent pattern to the red points.

Since I’m sampling uniformly throughout the square, there’s no reason to plot both where the identity holds and where it doesn’t. So for my next plot I just plotted where the identity fails.

Red stars scattered randomly among green dots

That’s a little  clearer. To make it clearer, I rerun the plot with 10,000 samples. (That’s the total number of random samples tested, about 16% of which are plotted.)

Red stars scattered randomly among green dots

Then to improve the resolution even more, I increased the number of samples to 1,000,000.

Red stars scattered randomly among green dots

It might not be too hard to work out analytically the regions where the identity holds and doesn’t hold. The blue regions above look sorta like parabolas, which makes sense because square roots are involved. And these parabola-like curves has ziggurat-like embellishments, which makes sense because floors are involved.

Update

Define

f(z) = ⌊√⌊z⌋ ⌋ − ⌊√z⌋.

The graphs above plotted the set of zs for which f(z) ≠ 0. We can see more information by plotting f itself. The function only takes on five possible values: 0, 1, i, −1, and −i. I plotted these as white, orange, green, blue, and red.

Related posts

[1] You could extend the meaning of “integer” to a complex number with integer real and imaginary parts, commonly known as a Gaussian integer.

[2] I suspect Mathematica’s definition is not common. I don’t know of any other source that even defines the floor of a complex number, much less defines it as Mathematica does.

[3] You may be wondering what’s going on with the square root. Complex numbers have two square roots, so which one(s) are we using? I’m doing what Python does, taking the principal branch. That is, using the analytic continuation of square root that extends from the real axis to the complex plane with the negative real axis excluded.