polar coordinates – Calculus VII

Pi and Python: how 22/7 morphs into 355/113

This is a brief return to the topic of Irrational Sunflowers. The sunflower associated with a real number ${a}$ is the set of points with polar coordinates ${r=\sqrt{k}}$ and ${\theta = a(2\pi k)}$ , ${k=1, 2, \dots, n}$ . A sunflower reduces to ${n}$ equally spaced rays if and only if ${a}$ is a rational number written in lowest terms as ${m/n}$ .

Here is the sunflower of ${a=\pi}$ of size ${n = 10000}$ .

Seven rays emanate from the center because ${\pi \approx 22/7}$ , then they become spirals, and spirals rearrange themselves into 113 rays because ${\pi \approx 355/113}$ . Counting these rays is boring, so here is a way to do this automatically with Python (+NumPy as np):

a = np.pi
n = 5000
x = np.mod(a*np.arange(n, 2*n), 1)
np.sum(np.diff(np.sort(x)) > 1/n)

This code computes the polar angles of sunflower points indexed ${5000\le k<10000}$ , sorts them and counts the relatively large gaps between the sorted values. These correspond to the gaps between sunflower rays, except that one of the gaps gets lost when the circle is cut and straightened onto the interval ${[0, 2\pi)}$ . So the program output (112) means there are 113 rays.

Here is the same sunflower with the points alternatively colored red and blue.

The colors blur into purple when the rational approximation pattern is strong. But they are clearly seen in the transitional period from 22/7 approximation to 355/113.

How many points would we need to see the next rational approximation after 355/113?
What will that approximation be? Yes, 22/7 and 355/113 and among the convergent of the continued fraction of ${\pi}$ . But so is 333/106 which I do not see in the sunflower. Are some convergents better than others?

Finally, the code I used to plot sunflowers.

import numpy as np
import matplotlib.pyplot as plt
a = np.pi
k = np.arange(10000)
r = np.sqrt(k)
t = a*2*np.pi*k 
plt.axes().set_aspect('equal')
plt.plot(r*np.cos(t), r*np.sin(t), '.')
plt.show()

Irrational sunflowers

A neat way to visualize a real number ${\alpha}$ is to make a sunflower out of it. This is an arrangement of points with polar angles ${2 \pi \alpha k}$ and polar radii ${\sqrt{k}}$ (so that the concentric disks around the origin get the number of points proportional to their area). The prototypical sunflower has ${\alpha=(\sqrt{5}+1)/2}$ , the golden ratio. This is about the most uniform arrangement of points within a disk that one can get.

But nothing stops us from using other numbers. The square root of 5 is not nearly as uniform, forming distinct spirals.

The number ${e}$ begins with spirals, but quickly turns into something more uniform.

The number ${\pi}$ has stronger spirals: seven of them, due to ${\pi\approx 22/7}$ approximation.

Of course, if ${\pi}$ was actually ${22/7}$ , the arrangement would have rays instead of spirals:

What if we used more points? The previous pictures have 500 points; here is ${\pi}$ with ${3000}$ . The new pattern has 113 rays: ${\pi\approx 355/113}$ .

Apéry’s constant, after beginning with five spirals, refuses to form rays or spirals again even with 3000 points.

Apéry's constant, 3000 points — Apéry’s constant, 3000 points

The images were made with Scilab as follows, with an offset by 1/2 in the polar radius to prevent the center from sticking out too much.

n = 500
alpha = (sqrt(5)+1)/2
r = sqrt([1:n]-1/2)  
theta = 2*%pi*alpha*[1:n]
plot(r.*cos(theta), r.*sin(theta), '*');
set(gca(), "isoview", "on")

Winding map and local injectivity

The winding map ${W}$ is a humble example that is conjectured to be extremal in a long-standing open problem. Its planar version is defined in polar coordinates ${(r,\theta)}$ by

$\displaystyle (r,\theta) \mapsto (r,2\theta)$

All this map does it stretch every circle around the origin by the factor of two — tangentially, without changing its radius. As a result, the circle winds around itself twice. The map is not injective in any neighborhood of the origin ${r=0}$ .

The 3D version of the winding map has the same formula, but in cylindrical coordinates. It winds the space around the ${z}$ -axis, like this:

In the tangential direction the space is stretched by the factor of ${2}$ ; the radial coordinate is unchanged. More precisely: the singular values of the derivative matrix ${DW}$ (which exists everywhere except when ${r=0}$ ) are ${2,1,1}$ . Hence, the Jacobian determinant ${\det DW}$ is ${2}$ , which makes sense since the map covers the space by itself, twice.

In general, when the singular values of the matrix ${A}$ are ${\sigma_1\ge \dots \ge\sigma_n}$ , the ratio ${\sigma_n^{-n} \det A}$ is called the inner distortion of ${A}$ . The word “inner” refers to the fact that ${\sigma_n}$ is the radius of the ball inscribed into the image of unit ball under ${A}$ ; so, the inner distortion compares this inner radius of the image of unit ball to its volume.

For a map, like ${W}$ above, the inner distortion is the (essential) supremum of the inner distortion of its derivative matrices over its domain. So, the inner distortion of ${W}$ is ${2}$ , in every dimension. Another example: the linear map ${(x,y)\mapsto (3x,-2y)}$ has inner distortion ${3/2}$ .

It is known that there is a constant ${K>1}$ such that if the inner distortion of a map ${F}$ is less than ${K}$ almost everywhere, the map is locally injective: every point has a neighborhood in which ${F}$ is injective. (Technical part: the map must be locally in the Sobolev class ${W^{1,n}}$ .) This was proved by Martio, Rickman, and Väisälä in 1971. They conjectured that ${K=2}$ is optimal: that is, the winding map has the least inner distortion among all maps that are not locally injective.

But at present, there is still no explicit nontrivial lower estimate for ${K}$ , for example we don’t know if inner distortion less than ${1.001}$ implies local injectivity.

Convexity of polar curves

Everybody knows the second derivative test for the convexity of Cartesian curves ${y=y(x)}$ . What is the convexity test for polar curves ${r=r(\theta)}$ ? Google search brought up Robert Israel’s answer on Math.SE: the relevant inequality is

$\displaystyle \mathcal{C}[r]:=r^2+2(r')^2-r\,r''\ge 0 \ \ \ \ \ (1)$

But when using it, one should be aware of the singularity at the origin. For example, ${r=1+\cos \theta}$ satisfies $\mathcal{C}[r] = 3(1+\cos \theta)\ge 0$ but the curve is not convex: it’s the cardioid.

The formula (1) was derived for ${r> 0}$ ; the points with ${r=0}$ must be investigated directly. Actually, it is easy to see that when ${r }$ has a strict local minimum with value ${0}$ , the polar curve has an inward cusp and therefore is not convex.

As usual, theoretical material is followed by an exercise.

Exercise: find all real numbers ${p}$ such that the polar curve ${r=(1+\cos 2\theta)^p}$ is convex.

All values ${p>0}$ are ruled out by the cusp formed at ${\theta=\pi/2}$ . For ${p=0}$ we get a circle, obviously convex. When ${p<0}$ , some calculations are in order:

$\displaystyle \mathcal{C}[r] = (1+4p+4p^2+(1-4p^2)\cos 2\theta)(1+\cos 2\theta)^{2p-1}$

For this to be nonnegative for all ${\theta}$ , we need ${1+4p+4p^2\ge |1-4p^2|}$ . Which is equivalent to two inequalities: ${p(4+8p) \ge 0}$ and ${2+4p\ge 0}$ . Since ${p<0}$ , the first inequality says that ${p\le -1/2}$ , which is the opposite of the second.

Answer: ${p=0}$ and ${p=-1/2}$ .

This exercise is relevant to the problem from the previous post Peanut allergy and Polar interpolation about the search for the “right” interpolation method in polar coordinates. Given a set of ${(r,\theta)}$ values, I interpolated ${(r^{1/p},\theta)}$ with a trigonometric polynomial, and then raised that polynomial to power ${p}$ .

If the given points had Cartesian coordinates ${(\pm a,0), (0,\pm b)}$ then this interpolation yields ${r=(\alpha+\beta \cos \theta)^p}$ where ${\alpha,\beta}$ depend on ${a,b}$ and satisfy ${\alpha>|\beta|}$ . Using the exercise above, one can deduce that ${p=-1/2}$ is the only nonzero power for which the interpolated curve is convex for any given points of the form ${(\pm a,0), (0,\pm b)}$ .

In general, curves of the form ${r=P(\theta)^{-1/2}}$ , with ${P}$ a trigonometric polynomial, need not be convex. But even then they look more natural than their counterparts with powers ${p\ne -1/2}$ . Perhaps this is not surprising: the equation ${r^2P(\theta)=1}$ has a decent chance of being algebraic when the degree of ${P}$ is low.

Random example at the end: ${r=(3-\sin \theta +2\cos \theta+\cos 2\theta-2\sin 2\theta)^{-1/2}}$

Peanut allergy and Polar interpolation

Draw a closed curve passing through the points ${(\pm 3,0)}$ and ${(0,\pm 1)}$ : more specifically, from (3,0) to (0,1) to (-3,0) to (0,-1) and back to (3,0).

I’ll wait.

… Done?

Okay, you may read further.

This is an interpolation problem. The geometry of the problem does not allow for curves of the form ${y=f(x)}$ . But polar graphs ${\rho=f(\theta)}$ should work. (Aside: in ${(x,y)}$ the variable usually taken to be independent is listed first; in ${(\rho,\theta)}$ it is listed second. Is consistent notation too much to ask for? Yes, I know the answer…) Since ${f}$ must be ${2\pi}$ -periodic, it is natural to interpolate ${(3,0),(1,\pi/2),(3,\pi),(1,3\pi/2)}$ by a trigonometric polynomial. The polynomial ${f(\theta)=2+\cos 2\theta}$ does the job. But does it do it well?

This is not what I would draw. My picture would be more like the ellipse with the given points as vertices:

The unnecessary concavity of the peanut graph comes from the second derivative ${f''}$ being too large at the points of minimum (and too small at the points of maximum). We need a function with flat minima and sharp maxima. Here is the comparison between ${f(\theta)=2+\cos 2\theta}$ (red) and the function that describes the ellipse in polar coordinates (blue):

I thought of pre-processing: apply some monotone function ${\varphi\colon (0,\infty)\rightarrow \mathbb R}$ to the ${\rho}$ -values, interpolate, and then compose the interpolant with ${\varphi^{-1}}$ . The function ${\varphi}$ should have a large derivative near ${0}$ , so that its inverse has a small derivative there, flattening the minima.

The first two candidates ${\varphi(r)=\log r}$ and ${\varphi(r)=1/r}$ did not improve the picture enough. But when I tried ${\varphi(r)=1/r^2}$ , the effect surprised me. Interpolation between ${(1/9,0),(1,\pi/2),(1/9,\pi),(1,3\pi/2)}$ yields ${(5-4\cos 2\theta)/9}$ , which after application of ${\varphi^{-1}}$ produces ${\rho = 3/\sqrt{5-4\cos 2\theta}}$ — exactly the ellipse pictured above.

Next, I tried a less symmetric example: curve through ${(5,0);(0,2);(-3,0);(0,-2)}$ .

Again, the second approach produces a more natural looking curve.

Finally, a curve with five-fold symmetry, with ${\rho}$ ranging from ${2}$ to ${5}$ .

The interpolation and plotting were done in Scilab, by invoking the functions given below with commands polar1([3,1,3,1]) or polar2([2,5,2,5,2,5,2,5,2,5]), etc. Because I considered equally spaced interpolation nodes, the coefficients of interpolated polynomials are nothing but the (inverse) discrete Fourier transform of the given values. First function interpolates the radius ${\rho}$ itself.

function polar1(R)
 clf()
 p = fft(R,1)  
 t = 0:.01:2*%pi
 rho = zeros(t)
 for i = 1:length(p)
     rho = rho + real(p(i)*exp(-%i*(i-1)*t))
 end
 polarplot(t,rho,style=5)
 t = resize_matrix(linspace(0,2*%pi,length(R)+1),1,length(R))
 plot(R.*cos(t), R.*sin(t), 'o')
endfunction

The second function includes preprocessing: it interpolates ${1/\rho^2}$ and then raises the interpolant to power ${-1/2}$ .

function polar2(R)
 clf()
 p = fft(R^(-2),1)
 t = 0:.01:2*%pi
 rho = zeros(t)
 for i = 1:length(p)
     rho = rho + real(p(i)*exp(-%i*(i-1)*t))
 end
 rho = max(rho,0.01*ones(rho))
 polarplot(t,rho^(-1/2),style=5)
 t = resize_matrix(linspace(0,2*%pi,length(R)+1),1,length(R))
 plot(R.*cos(t), R.*sin(t), 'o')
endfunction

Added: one more comparison, ${\rho^{1}}$ vs ${\rho^{-2}}$ vs ${\rho^{-3}}$ ; the last one loses convexity again.

And then to ${\rho^{-7}}$ :