Area of a triangle in the complex plane

I recently ran across an elegant equation for the area of a triangle in the complex plane with vertices at z1, z2, and z3. [1].

A(z_1, z_2, z_3) = \frac{i}{4} \, \left| \begin{matrix} z_1 & \bar{z}_1 & 1 \\ z_2 & \bar{z}_2 & 1 \\ z_3 & \bar{z}_3 & 1 \\ \end{matrix} \right|

This formula gives the signed area: the area is positive if the points are given in countclockwise order and negative otherwise.

I’ll illustrate the formula with a little Python code. Let’s generate a random triangle.

    import numpy as np

    np.random.seed(20221204)
    r = 100*np.random.random(6)
    z1 = r[0] + 1j*r[1]
    z2 = r[2] + 1j*r[3]
    z3 = r[4] + 1j*r[5]

Here’s what our triangle looks like plotted.

Now let’s calculate the area using the formula above and using Heron’s formula.

    
    def area_det(z1, z2, z3):
        det = 0
        det += z2*z3.conjugate() - z3*z2.conjugate()
        det -= z1*z3.conjugate() - z3*z1.conjugate()
        det += z1*z2.conjugate() - z2*z1.conjugate()
        return 0.25j*det
    
    def area_heron(z1, z2, z3):
        a = abs(z1-z2)
        b = abs(z2-z3)
        c = abs(z3-z1)
        s = 0.5*(a + b + c)
        return np.sqrt(s*(s-a)*(s-b)*(s-c))
        
    print(area_heron(z1, z2, z3))
    print(area_det(z1, z2, z3))

This prints -209.728 and 209.728. The determinate gives a negative area because it was given the points in clockwise order.

[1] Philip J. Davis. Triangle Formulas in the Complex Plane. Mathematics of Computation. January 1964.

Conformal map between square and disk

Conformal maps transform one region into another while preserving angles. You might solve a PDE, for example, by mapping it to a standard region, solving it there, then mapping the solution back to the original region.

Some tasks are easier to do in a square and others in a disk, so it’s clearly useful to be able to conformally map between squares and disks. The Riemann mapping theorem tells us this can be done, but it doesn’t tell us how. Two gentlemen figured out how to map between squares (and more general polygons) and disks in the 1860s: Hermann Schwarz and Elwin Christoffel. Schwarz is known for many different results in analysis, including the topic of the previous post, the conformal map from an ellipse to the unit disk. Christoffel is best known for Christoffel symbols, building blocks of tensors.

Here’s a plot showing how the Schwarz-Christoffel transformation from the square [-1, 1] × [-1, 1] to the unit disk transforms Cartesian grid lines.

Here’s another plot, this one showing how the grid lines for polar coordinates on the disk pull back to curves on the square.

Equations

The equation for the function from the square to the disk is

f(z) = \frac{1-i}{2} \, \text{sd}\left(\frac{1+i}{2} K z;\frac{1}{2}\right)

where sd is a Jacobi elliptic function with parameter 1/2 [1]. The constant K is the complete elliptic function of the first kind, evaluated at 1/2. In symbols, K = K(1/2).

The inverse function has equation

g(w) = (1-i) - \frac{1-i}{K} F\left(\cos^{-1}\left(\frac{1-i}{\sqrt{2}} w\right);\frac{1}{2} \right )

Here F is the incomplete elliptic function of the first kind. For more background, see this post on kinds of elliptic integrals.

Peirce’s projection

Charles Sanders Peirce used the conformal map of the disk to the square to create the “Peirce quincuncial projection” map. This is a conformal (i.e. angle-preserving) map that represents the globe on a square. The diamond shape in the middle is the image of the equator. The mapping is singular at the south pole.

Charles Sanders Peirce's quincuncial project

Peirce named the map after the quincunx pattern of the poles. This obscure word refers to the pattern of dots on the five face of a standard six-sided die.

⚄

Related posts

[1] “Parameter” is being used in a technical sense here. There are two conventions for describing the parameterization of elliptic functions and elliptic functions and elliptic integrals, and here we are using the parameter called The Parameter, commonly denoted m. There’s another convention that uses the elliptic modulus k, and the connection between them is that m = k².

Conformal map of ellipse interior to a disk

This post will present the conformal map between the interior of an ellipse and the unit disk.

Given an ellipse centered at the origin with semi-major axis a and semi-minor axis b. Will will assume without loss of generality that a² – b² = 1 and so the foci are at ±1.

Hermann Schwarz published the conformal map from the ellipse to the unit disk in 1869 [1, 2].

The map is given by

f(z) = \sqrt{k}\,\, \text{sn} \left( \frac{2K}{\pi} \sin^{-1}(z) \right)

where sn is the Jacobi elliptic function with parameter k². The constants k and K are given by

\begin{align*} q &= (a + b)^{-4} \\ k &= \left(\frac{\theta_2(q)}{\theta_3(q)}\right)^2 \\ K &= \frac{\pi}{2} \theta_3(q)^2 \end{align*}

where θ2 and θ3 are theta constants, the value so the theta functions θ2(z, q) and θ3(z, q) at z = 1.

Conformal maps to the unit disk are unique up to rotation. The map above is the unique conformal map preserving orientation:

\begin{align*} f(0) &= 0 \\ f(\pm a) &= \pm 1 \\ f(\pm bi) &= \pm i \\ \end{align*}

Inverse map

The inverse of this map is given by

g(w) = \sin\left(\frac{\pi}{2K} \,\, \text{sn}^{-1} \left(\frac{w}{\sqrt{k}}\right) \right)

The inverse of the sn function with parameter m can be written in terms of elliptic integrals.

\text{sn}^{-1}(z; m) = F\left(\sin^{-1}(u; m)\right)

where F is the incomplete elliptic integral of the first kind and m is the parameter of sn and the parameter of F.

Plot

I wanted to illustrate the conformal map using an ellipse with aspect ratio 1/2. To satisfy a² – b² = 1, I set a = 2/√3 and b = 1/√3. The plot at the top of the post was made using Mathematica.

Related posts

[1] H. A. Schwarz, Über eigige Abbildungsaufgaben, Journal für di reine und angew. Matheamatik, vol 70 (1869), pp 105–120

[2] Gabor Szegö. Conformal Mapping of the Interior of an Ellipse onto a Circle. The American Mathematical Monthly, 1950, Vol. 57, No. 7, pp. 474–478

Eliminating a Bessel function branch cut

In an earlier post I looked in detail at a series for inverse cosine centered at 1. The function arccos(z) is multivalued in a neighborhood of 1, but the function

arccos(z) / √(2 – 2z)

is analytic in a neighborhood of 1. We cancel out the bad behavior of inverse cosine at 1 by dividing by another function with matching bad behavior.

We can do something similar with Bessel functions. In general, the Bessel function Jν has a branch cut from -∞ to 0. But the function

Jν / zν

is entire, i.e. analytic everywhere in the complex plane. From a certain perspective—namely the perspective advocated by Bille Carlson (1924–2013)—we got the definition of Bessel functions (and other functions) wrong. Carlson’s book Special Functions of Applied Mathematics shows how many special functions can be defined in terms of a smaller number of functions with fewer branch cuts. This simplification revealed symmetries that hadn’t been noticed before.

Here’s a plot of J1/3(z). Height represents magnitude and color represents phase.

This was produced using the Mathematica command

    ComplexPlot3D[BesselJ[1/3, z], {z, -1 - I, 1 + I}]

The abrupt change from purple to yellow along the negative real axis shows the discontinuity in the phase at the branch cut. You can also see a hint of the cusp at the origin.

Here’s a plot of J1/3(z) / z1/3.

This was made using

    ComplexPlot3D[BesselJ[1/3, z]/z^(1/3), {z, -1 - I, 1 + I}]

The white streak is an artifact of the branch cuts for J1/3(z) and for z1/3. The branch cut in their ratio is unnecessary.

Related posts

Complex AGM

The arithmetic-geometric mean (AGM) of two non-negative real numbers a and b is defined as the limit of the iteration starting with a0 = a and b0 = b and

an+1 = ½ (an + bn)
bn+1 = √(an bn)

for n > 0. This sequence converges very quickly and is useful in numerical algorithms. The limit can be expressed in terms of an elliptic function, and that elliptic function can then be related to other functions. See, for example, this post for how the AGM can be used to compute logarithms to arbitrary precision.

Since the AGM is useful in computing special functions, and we’re often interested in evaluating special functions at complex values, it’s natural to want to evaluate the AGM for complex numbers.

But we immediately run into a difficulty: which square root do we pick to find the new b?

For a non-negative real number x, √x is defined to be the non-negative real number y such that y² = x. But for more general values of x we have to choose a branch of the square root function. Which branch should we use in the AGM?

Often when we need to take the square root of complex numbers we can use the “principal branch,” the branch that gives positive values for positive real inputs and extends to the rest of the complex plane with the negative axis removed. If you compute the square root of a complex number in software, as we’ll do below, this is probably the value you’ll get by default.

But it turns out we cannot simply pick the principal branch. Gauss discovered two centuries ago that the right branch to take could vary at each step [0]. What does that even mean? How do we know we’ve made the “right” choice?

For one thing, we want our iteration to converge. And for another, we’d like it to converge to something non-zero if we start with non-zero inputs [1]. The right choice will guarantee this [2].

So what is this right choice? We provisionally update a and b as above, using either square root for b, and keep the value of b if

|ab| ≤ |a + b|

or

|ab| = |a + b|

and the imaginary part of b/a is positive. In other words, chose the possible value of b that’s closer to a, and use the imaginary part of b/a as a tie breaker.

Here is the AGM implemented in Python, using the right square root at each step.

    def right(a, b):
        d = abs(a + b) - abs(a - b)
        return d > 0 or d == 0 and (b/a).imag > 0

    def agm(a, b):
        while abs(a-b) > 1e-14:
            a1 = 0.5*(a+b)
            b1 = np.sqrt(a*b)
            if not right(a1, b1):
                b1 = -b1
            a, b = a1, b1
        return a1

The test d == 0 should make you concerned. Testing floating point numbers for exact equality with zero is seldom the right thing to do, but we’ll gloss over numerical issues for this post.

Here’s an example, giving the a values of the iteration starting with 7 + 30i and 20 + 22i. The iteration converges in four steps.

    13.500000000000000 + 26.000000000000000j
    13.784944719026262 + 26.397404494892115j
    13.783557503026870 + 26.395953326186888j
    13.783557473769877 + 26.395953309190112j

In this example, the right root is the principal root every time. To find an example where the other root is chosen, I generated random starting points until I found one that took an alternate root.

Here are the values of b starting with -1.654 – 1.178i and 2.244 – 1.956i. An asterisk indicates that the principal root was not the right root.

     0.2328790598285062 - 0.728412421988127j
    -0.6829999999999998 - 0.589000000000000j
    -0.2254063569280081 - 0.799311791579126j *
    -0.2250604700857468 - 0.658706210994063j
    -0.2261796153095098 - 0.725905503054624j *
    -0.2252334135068774 - 0.729009001286595j
    -0.2257078598391289 - 0.727456168426402j *
    -0.2257065144081936 - 0.727457252170610j
    -0.2257071871240875 - 0.727456710298747j *
    -0.2257071871236612 - 0.727456710298506j
    -0.2257071871238743 - 0.727456710298626j *
    -0.2257071871238744 - 0.727456710298626j

More AGM posts

[0] David Cox. The Arithmetic-Geometric Mean of Gauss. L’Énseignement Mathématique, 30 (1984), p. 275–330.

[1] If a = –b we get zero immediately. But if a and b are both not zero, and a does not equal –b, then taking the right square root at each iteration gives us what we want.

[2] Interestingly, you could make a finite number of wrong choices and still end up with something that might converge, albeit to a different value. This gives you a different branch of the AGM.

Keyhole contour integrals

The big idea

The Cauchy integral theorem says that the integral of a function around a closed path in the complex plane depends only on the poles of the integrand inside the path. You can change the path itself however you like as long as you don’t change which poles are inside. This observation is often used to compute real integrals using complex analysis.

Suppose you want to integrate a function f along (some portion of) the real line, and f extends to a function in the complex plane that is analytic except at poles. You may be able to evaluate your real integral using a complex contour, or more commonly, a limit of contours. You show that in some limit, the contribution to the integral along the parts of the contour you don’t need goes to zero, and the rest of the contour approaches the part you wanted to integrate over initially.

Keyhole contours

There area a handful of integration contours that come up most frequently, and the keyhole is one of them. The initial motivation for this post was the cover of the fourth edition of A Course of Modern Analysis by Whittaker and Watson. I’ve written a couple posts lately about the story behind images on book covers, and this is another post in that series.

E. T. Whittaker & G. N. Watson, A Course of Modern Analysis

The most recent edition has a more prosaic cover with only words and no image. However, the latest addition is much more attractive inside, having been rewritten in LaTeX.

The image on the cover is a keyhole contour. Typically the slot in the contour runs along the real axis, but the cover rotated the image 45 degrees to make it more visually appealing. The contour is used on page 118 of the fourth edition to integrate rational functions along the positive real axis [1].

Contour integration drawing

If the function you’re integrating has poles only at the locations marked by stars, then you can evaluate the integral around the contour by computing the residues of the integrand at these points. Now suppose you change the contour by making the outer circle larger, letting its radius R go to infinity, and making the inner circle shrink, letting its radius ε go to zero. If the integral along these two circular segments goes to zero in the limit, you’re left with the integral along the positive real axis.

Hankel

The keyhole contour is sometimes called the Hankel contour because Hermann Hankel used it in the 19th century to investigate special functions such as the gamma function and the eponymous Hankel functions. The first post in this series of book cover commentaries mentioned a book that has a plot of a Hankel function on the cover.

Related posts

[1] Specifically, Whittaker and Watson show how to compute the integral of xa-1 Q(x) from 0 to infinity, where Q is a rational function with no positive real zeros, and the limit of xa Q(x) is zero both as x goes to zero and as x goes to ∞.

NASA and conformal maps

A couple years ago I wrote about how NASA was interested in regions bounded by curves of the form

\left| \frac{x}{A} \right|^\alpha + \left| \frac{y}{B} \right|^\beta = 1

For example, here’s a plot for A = 2, B = 1, α = 2.5 and β = 6.

That post mentioned a technical report from NASA that explains why these shapes are important in application, say for beams and bulkheads, and gives an algorithm for finding conformal maps of the unit disk to these shapes. These shapes are also used in graphic design, such as squircle buttons on iPhones.

Note that these shapes are not rounded rectangles in the sense of a rectangles modified only at the corners. No segment of the sides is perfectly flat. The curvature smoothly decreases as you move away from the corners rather than abruptly jumping to 0 as it does in a rectangle rounded only at the corners. Maybe you could call these shapes continuously rounded rectangles.

Conformal maps

Conformal maps are important because they can map simple regions to more complicated regions in way that locally preserves angles. NASA might want to solve a partial differential equation on a shape such as the one above, say to calculate the stress in a cross section of a beam, and use conformal mapping to transfer the equation to a disk where the calculations are easier.

Coefficients

The NASA report includes a Fortran program that computes the coefficients for a power series representation of the conformal map. All the even coefficients are zero by reasons of symmetry.

Coefficients are reported to four decimal places. The closer the image is to a circle, i.e. the closer α and β are to 2 and the closer A is to B, the fewer non-zero coefficients the transformation has. You could think of the number of non-zero coefficients as a measure how hard the transformation has to work to transform the disk into the desired region.

There are interesting patterns in the coefficients that I would not have noticed if I had not had to type the coefficients into a script for the example below. Maybe some day I’ll look more into that.

Conformal map example

The following example uses A = 2, B = 1, and α = β = 5.

We start with the following disk with polar coordinate lines drawn.

And here is the image of the disk under the conformal map giving by the coefficients in the NASA report.

Note that the inner most circles from the domain are mapped into near circles, but as you move outward the circles are distorted more, approaching the shape of the boundary of the range.This is because in the limit, analytic functions map disks to disks. So a tiny circle around the origin will always map to an approximate circle, regardless the target region.

Note also that although the circles and lines of the domain are warped in the range, they warp in such a way that they still intersect at right angles. This is the essence of a conformal map.

Inverse map

The NASA report only gives a way to compute the transformation from the disk to the target region. Typically you need to be able to move both ways, disk to target and target to disk. However, since the transformation is given by a power series, computing the inverse transformation is a problem that has been solved in general. This post gives an explanation and an example.

Related posts

Transformations of Olympic rings

The previous post gave the details of how Möbius transformations

m(z) = (az + b)/(cz + d)

transform circles. The image of a circle under a Möbius transformation is either a line or a circle, and in our examples the image will always be a circle.

We start with an approximation of the Olympic rings

 

 

and then apply various transformations.

First we apply (z + 9)/(z + 10):

 

This transformation doesn’t distort the rings too much because the rings are far from the singularity at z = -10. The blue ring is distorted the most because it’s closest to the singularity.

Next we apply (iz + 7)/(z + 8).  Now we’ve moved the singularity a little bit closer, and we caused some rotation by using i as the coefficient of z in the numerator.

 

 

Finally, let’s pick a transformation with its singularity inside the green circle, (z + 4)/(z -3 + i).

 

 

Now we see a lot of distortion. The blue and gold rings that were outside the green ring are brought inside. The complex plane is turned inside-out because of the singularity.

As before, the black and red rings still intersect the green ring. But what you can’t tell from the graph is that the portion of these rings that was inside has been turned outside and vice versa.

Related posts

Circles and lines under a Möbius transformation

This post will revisit a previous post in more detail. I’ve written before about how Möbius transformations map circles and lines to circles and lines. In this post I’d like to explore how you’d calculate which circle or line a given circle or line goes to. Given an equation of a line or circle, what is the equation of its image?

Circles to circles

A while back I wrote a post called Circles to Circles, demonstrating how Möbius transformations map circles and lines to circles and lines. If you think of a line as a degenerate circle, a circle of infinite radius, then you could simply say Möbius transformations map circles to circles.

For example I showed in that post that the transformation

m(z) = (z + 2)/(3z + 4)

maps these circles

to these circles

and that these lines

go to these circles

Circle or line image

Given the equation of a circle or a line, how would you know whether its image is a circle or a line?

A Möbius transformation

m(z) = (az + b) / (cz + d)

has a singularity at z = –d/c. If a line or curve goes through the point –d/c then its image is a line. Otherwise the image is a circle.

If the line or curve does go through –d/c, so that the image is a line, it’s easy to find an equation of the line. Pick two points, z1 and z2 on the original curve or line. Let w1m(z1) and w2m(z2). Then find the equation of a line through w1 and w2. For example, a parametric form is

w1 + t(w2w1)

for all real t.

Equation of the image circle

Now suppose we start with a circle centered at γ with radius r. We apply

m(z) = (az + b) / (cz + d)

as above, but we assume our circle does not pass through –d/c since that case is covered above. With this assumption we can conclude that the image of our circle is another circle, and the remaining task is to find an equation of that circle.

We can decompose our transformation m into several simpler steps—translation, reciprocal, rotation, and dilation—and see what happens to the center and radius under each step. The hardest step is taking the reciprocal, and this is covered in the previous post. That post was a lemma for this post.

We can apply m to z in four steps:

  1. Add d/c.
  2. Take the reciprocal.
  3. Multiply by e = (bcad)/c².
  4. Add a/c.

We know it’s OK to take the reciprocal in the second step because our circle does not contain the point z = –d/c.

We can break step 3 into (3a) multiplying by |e| and (3b) multiplying by e/|e|. Step 3a is a dilation and step 3b is a rotation.

Here’s a Python program using the steps above to find the center and radius of a transformed circle.

    def transform(a, b, c, d, center, radius):
        # step 1
        center += d/c

        # step 2
        temp = center
        center = np.conj(center) / (abs(center)**2 - radius**2)
        radius /= abs(abs(temp)**2 - radius**2)

        # step 3
        e = (b*c - a*d)/c**2
        center *= e
        radius *= abs(e)

        # step 4
        center += a/c

        return (center, radius)

Computing zeta at even numbers

Last year I wrote several posts about computing ζ(3) where ζ is the Riemann zeta function. For example, this post.

It happens that ζ can be evaluated in closed form at positive even arguments, but there’s still a lot of mystery about zeta at positive odd arguments.

There’s a way to derive ζ(2n) using contour integration. But you don’t need to be clever about constructing contours if you start from the following result which is derived using contour integration.

Suppose f is analytic except for poles at zk and that the sum

\sum_{n \ne z_k} f(n)

converges, where the sum is over all integer n except poles of f. Then the sum can be computed from residues at the poles:

 \sum_{n \ne z_k} f(n) = -\sum_{k} \text{Res}\left( f(z) \pi \cot \pi z; z_k \right)

Here Res stands for “residue.” The residue of a function g(z) at z0 is defined as the coefficient c-1 of (zz0)-1 in the Laurent series for g. Equivalently,

 c_{-1} = \lim_{z\to z_0} (z - z_0) g(z)

This means we can evaluate infinite sums by taking a limit that may be a simple application of L’Hôpital’s rule.

The theorem above could be used to calculate a lot more than values of the Riemann zeta function by choosing various functions f, but for our purposes f(z) will be z-2m for some positive integer m.

\zeta(2m) = \sum_{n=1}^\infty \frac{1}{n^{2m}} = \frac{1}{2} \sum_{\stackrel{n=-\infty}{n\ne 0} }^\infty \frac{1}{n^{2m}} = -\frac{1}{2}\,\text{Res}\left( \frac{\pi \cot \pi z}{z^{2m}}; 0 \right )

Why can’t this same trick be used for evaluating the Riemann zeta function at odd integer arguments? The second equality above fails. For even powers of n, the sum over the positive integers is half the sum over all integers. For odd powers, the sum over all integers is zero.

If we write out the Laurent series for π cot πz then we can read off the residues of π cot πz/z2m.

\pi \cot \pi z = \frac{1}{z}-\frac{\pi ^2 z}{3}-\frac{\pi ^4 z^3}{45}-\frac{2 \pi ^6 z^5}{945} - \cdots

When we divide by z² the coefficient of z becomes the coefficient of z-1, i.e. the residue. This residue is -π²/3, and so negative one half of this is π²/6, i.e. ζ(2) = π²/6.

The same line of reasoning shows that ζ(4) = π4/90 and ζ(6) = π6/945.

Finding the values of ζ(2m) in general is no harder than finding the power series for cotangent, which is not trivial but not insurmountable either. See this post.

The final result is that

\zeta(2m) = (-1)^{m-1} 2^{2m-1} \pi^{2m} \frac{B_{2m}}{(2m)!}