Phyllotaxis and Fibonacci

Epicycles

Any orbit can be represented given enough epicycles.

The Geocentric Model

The geocentric model was wonderful for human self aggrandizement.  We belonged at the centre of the universe and everything revolved around us.  That early humans believed this isn't that surprising.  After all the Sun, moon, and stars do seem to move in circles around us (albeit different ones).  What's more interesting is how we attempted to cling on to this theory in the light of a) conflicting evidence, and b) a far better explanation.

David Deutsch in The Fabric of Reality puts forward the view that the point of a scientific theory is to explain.  The more a theory explains - whilst still remaining consistent with observable facts - the better it is.  This is a philosophical justification for Occam's Razor.  Why is an explanation with fewer postulates better than one with more postulates? because it leaves less unexplained!

The earliest theory for why the planets did not move in simple circular orbits like the Sun, moon, and stars, was probably that they were sentient beings with their own free will.  The word planet is derived from the ancient Greek for wanderer, and the planets known in antiquity share names with classical gods.  But when formulae were found for their paths (based on epicycles, which I'll come to later) the sentient wanderer theory had to bite the dust.  In light of the new facts, the sentient wanderer theory became "the planets don't move like the Sun and the moon because they are sentient beings that choose their own paths, but they always choose these ones".  This leaves more unexplained than just saying "the planets paths are given by these formulae" because as well as failing to explain why these formulae were chosen it also fails to explain what the gods are doing up there.  In terms of Occam, the additional postulate of there being gods up there adds nothing and so should go.

Epicycles


In the epicycle theory the paths were not circular orbits, but circular orbits around circular orbits.  Later, when more accurate measurements were taken, it became necessary to model them by circular orbits around circular orbits around circular orbits.  If we simplify things slightly by taking the planets to all lie in a plane then we can use the complex numbers to represent positions.  A circular orbit taking time $T$ around $a_0$ becomes
$$
a_0 + a_1e^{i\omega t}
$$
where $a_1$ is a complex number defining the radius and phase of the orbit and $\omega = 2\pi/T$.  For circular orbits around circular orbits lets restrict ourselves to double the angular frequency.  The equation becomes
$$
a_0 + a_1e^{i\omega t} + a_2e^{2 i \omega t}
$$
And for circular orbits about circular orbits around circular orbits if we restrict the frequency to three times $\omega$ we get
$$
a_0 + a_1e^{i\omega t} + a_2e^{2 i \omega t} + a_3e^{3 i \omega t}
$$
You can see where this is going, if we do not limit ourselves to a finite number of epicycles the trajectory becomes $\sum_{n=0}^\infty a_ne^{i  \omega n t}$.  These belong to the vector space of continuous and cyclic functions $[0,T] \to \mathbb{C}$.  (Cyclic just means $f(0)=f(T)$.)

But astronomers only make a finite number of measurements per orbit.  So let's restrict the vector space to  functions $[t_0, ..., t_{N-1}] \to \mathbb{C}$ (where $t_n = \frac{n}{N}T$).  This has a simple inner product
$$
\left<f, g\right> = \frac{1}{N}\sum_{n=0}^{N-1}{f \left(t_n\right)^*g\left(t_n\right)}
$$
And it's easy to check that $\{e^{i  \omega n t} : 0 \le n \lt N\}$ forms an orthonormal basis under this inner product.  What does this mean?  It means that any periodic orbit around us, however squiggly, can be represented by N epicycles leading to the formula below.  And if you only make N measurements the model predicts the measurements perfectly!
$$
\begin{align}
f(t) &=
a_0 + a_1e^{i\omega t} + a_2e^{2 i \omega t} + ... + a_{N-1}e^{(N-1) i \omega t} \\
\text{where,} & \\
a_n &= \left<e^{i \omega n t}, f\left(t\right)\right>
\end{align}
$$

The Heliocentric Model

So it turns out that the theory of epicycles can explain any periodic trajectory.  But in general if you are going to make N observations you need N epicycles to explain them.  But this doesn't really explain anything, all it does is shift the problem.  Instead of needing to explain why we have these N measurements we are left needing to explain why we have these N epicycles.

At first this wasn't obvious because the first few epicycles provided better resolution than the measurements.  However, the Heliocentric theory motivated an effort to make more accurate measurements, which might distinguish between the two theories.  When these better measurements started supporting the Heliocentric model the establishment response wasn't to ditch the Geocentric model, but to add new epicycles!  It's easy to be unfair on the renaissance scientific establishment, looking back with the benefit of Fourier Analysis.  But even without the understanding that with enough epicycles you can generate any path to any desired accuracy, the principle of Occam's Razor should have been enough to choose between two theories.  One kept needing to grow as new facts came in, and the other remained the same size.

Code for the illustration


#!/usr/bin/env python2
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

def inner_product(a,b):
    'Inner product of two complex N-vectors'
    return (a.conj() * b).sum() / len(b)

def get_nth_orthonormal_vector(n, N):
    'Return the nth orthornomal N-vector e^{2 \pi i n / N}'
    return np.exp(2 * np.pi * n / N * 1j * np.arange(N))

def discrete_fourier_transform(a):
    'Return coefficients np.array([A_0,A_1,...]) relative to orthonormal basis'
    N = len(a)
    return np.array(map(lambda n: inner_product(get_nth_orthonormal_vector(n, N),a), range(N)))

def inverse_discrete_fourier_transform(A):
    'Return the original samples np.array([a_0,a_1,...]) given the orth. basis coefficients'
    N = len(A)
    return sum(map(lambda n: A[n] * get_nth_orthonormal_vector(n, N), range(N)))

# Create a silly square orbit z around 0 + j0
x = np.array(range(-10,10) + [10]*20 + range(10,-10,-1) + [-10]*20)
y = np.array([-10]*20 + range(-10,10) + [10]*20 + range(10,-10,-1))
z = x + 1j * y
N = len(z)

# Get the DFT of z
z_dft = discrete_fourier_transform(z)

# Verify iDFT of z_dft is z (allowing for rounding errors)
assert np.max(np.abs(inverse_discrete_fourier_transform(z_dft) - z)) < 1e-10

# Animate
fig = plt.figure()
ax = plt.axes(xlim=(-13,13), ylim=(-13,13))
ax.set_aspect('equal')
ax.axis('off')
components = map(lambda n: z_dft[n] * get_nth_orthonormal_vector(n, N), range(N))
circle_objs = map(lambda n: ax.plot([], [], color='b', alpha=0.5)[0], range(N))

def update(t):
    pos = 0j
    for n in [ -k/2 if (k&1) else k/2 for k in range(N) ]:
        circle = np.ndarray(N+1, dtype=complex)
        circle[:N] = components[n] + pos
        circle[N] = circle[0]
        circle_objs[n].set_data(circle.real, circle.imag)
        pos = circle[t]
    scat_obj = ax.scatter([pos.real], [pos.imag], color='b', marker='.', alpha=0.5)
    return tuple(circle_objs) + (scat_obj,)

ani = animation.FuncAnimation(fig, update, frames=N, interval=100, blit=True, repeat=True)
ani.save('output.gif', dpi=80, writer='imagemagick')
#plt.show()

Postscript

The algorithm used above is actually identical to that used in OFDM - which is the basis for all 4G/5G transmissions.  In OFDM it is used to convert between I/Q samples at baseband frequency and per-subcarrier I/Q coefficients.  Here it's used to convert between trajectory points and epicycle coefficients, but the maths is identical.

Comments