## Liebeck’s sequence: an applied mathematician writes

Our solution to Liebeck’s problem, and in particular our comments about how an applied mathematician might tackle it, provoked a response from… an applied mathematician, Prof. Nigel Mottram. Here, with permission, are his thoughts on the problem.

If I have to find a continuous analytic function which approximates this sequence of numbers (i.e. if I’m too lazy to work out if $n$ is closer to $3^m$ or $2\cdot3^m$) then this is what I’d do…

As an applied mathematician I know that I must
(a) not be afraid to experiment, and
(b) not be afraid to use previous results from any discipline if they can help my cause (suitably cited of course).

So… first I use Matlab to do some experiments (see the code at the end of this post) and I see that our sequence $f(n)$ is just a sequence of sequences that lie on lines with gradient 1 or 3. The sequence $f(n)$ changes its gradient if it hits the line $y=2x$ or $y=3x/2$.

This agrees with your result that $f(3^m)=2 \cdot f(3^m)$ and $f(2\cdot 3^m)=3^{m+1}$. I could have just used your result to get to this point of course…

So I now have a sequence of numbers which seem to oscillate around the line $y=7x/4$ (which is the average of $y=2x$ and $y=3x/2$)… so I’ll use this result to say that $f(n)=7n/4$ is my first approximation.

I’m an applied mathematician so I love Fourier and believe everything is essentially a sine wave. So now I try and think of a sinusoidal perturbation of this first order approximation. I realise that the amplitude of the perturbation is getting larger but it stays within (and attains) the bounds of $y=2x$ and $y=3x/2$ so the amplitude of my sine wave grows like $x/4$ (the maximum minus the average value, $2-7/4$). The period looks like it is growing as well. In fact I know that the sequence hits the line $y=2x$ whenever $x$ is of the form $3^m$… which means I need a sinusoidal function which hits its maximum whenever $\log(x)/\log(3)$ is an integer. I’m therefore looking at a better approximation which is

$f(n)= \displaystyle\frac{7n}{4} + \displaystyle\frac{n}{4}\cos\left(\displaystyle\frac{2\pi\log(n)}{\log(3)}\right)$

and gives us $f(100)=184$.

I could think about when the sequence reaches $y=3x/2$ instead and I get an approximation

$f(n)= \displaystyle\frac{7n}{4} - \displaystyle\frac{n}{4}\cos\left(\displaystyle\frac{2\pi(\log(n)-\log(2))}{\log(3)}\right)$

and this gives us $f(100)=198$.

The first of these approximations is better when I’m close to $y=2x$ and the latter when I’m close to $y=3x/2$… but if I don’t know this beforehand then maybe an average of these two gives a lower global error.

If we look at the first of these approximations we can rewrite it as

$f(n)= \displaystyle\frac{n}{4}\left( 7+\cos\left(\displaystyle\frac{2\pi\log(n)}{\log(3)}\right) \right).$

Now let’s change our coordinates slightly making a new coordinate $t=\log(n)$ and consider a rescaled sequence $F(n)=4f(n)/n$. This approximation is now

$F(n)=7+\cos\left(\displaystyle\frac{2\pi t}{\log(3)}\right)$

which certainly looks like the start of a Fourier series. If we take the original sequence and plot $F(n)=4f(n)/n$ against $\log(n)$ then we see a lovely periodic function which can be expressed as a Fourier series… it looks like we are home and dry.

But be careful…if we only take, let’s say, $k$ terms in our Fourier series (call the Fourier series $S_k$) there will be an error (say $\epsilon_k$) and our approximation will look like

$f(n)=\displaystyle\frac{n}{4}(S_k+\epsilon_k),$

and we see that the error in the reconstructed sequence will grow linearly with $n$.

So what we have found is that we can approximate the sequence $f(n)$ using an continuous analytic function but unfortunately the error of this approximation grows as n increases.

It would be far better to go back to the beginning and use the fact that the sequence is simply a sequence of linearly increasing numbers (with gradients 1 and 3) between the number $3^m$ and $2\cdot 3^m$.

(NJM)

The Matlab code used to generate these figures was:

clear all

a(1)=2;
n=1000;
for j=2:n
i=find(a(1:j-1)==j);
if isempty(i)==0;
a(j)=3*i;
else
a(j)=a(j-1)+1;
end
end

figure(1)
plot(1:n,a,’go’);
hold on
plot(1:n,1.5*(1:n),’r-‘)
plot(1:n,2*(1:n),’r-‘)
hold off

figure(2)
plot(1:n,a,’go’);
hold on
plot(1:n,1.5*(1:n),’r-‘)
plot(1:n,2*(1:n),’r-‘)
x=1:n/1000:n;
plot(x,1.75*x+0.25*x.*(cos(2*pi*(log(x)/log(3)))),’b-‘)
hold off

figure(3)
plot(1:n,a,’go’);
hold on
plot(1:n,1.5*(1:n),’r-‘)
plot(1:n,2*(1:n),’r-‘)
x=1:n/1000:n;
plot(x,1.75*x-0.25*x.*(cos(2*pi*((log(x)-log(2))/log(3)))),’b-‘)
hold off

figure(4)
t=log(1:n);
plot(t,4*a./exp(t),’go-‘);
hold on
plot(t,6,’r-‘)
plot(t,8,’r-‘)
hold off