I am closing the blog

I am closing the blog and I will be (slowly) moving the posts to another one which I have started on Github, written with Octopress.

Main reasons :

  • There are now advertisements inside the main body of the blogs on WordPress.com. How lame is that !?
  •   Octopress handles maths and code much better, you can write your post in markdown, its themes are fully customizable for free, and you don’t need an internet access to make your post drafts !
  • This blog can too easily be found when googling my name.

[Doesn’t work anymore, sorry, might fix that one day…] Customize your IPython notebook with CSS

In the coming IPython 1.0 (available on their github) you will be able to change the style of your notebook using a custom CSS file. Here are three attempts, which I will describe in details:

Creating a IPython profile

Credit and many thanks for this part go to Matthias Bussonier for his thorough explanatory blog on the subject.

Creating a profile is easy, you just type this line in a terminal:

ipython profile create yourprofilename

Note that you can send commands to the terminal from within your IPython notebook by writing %%bash in the first line of a code cell:

%%bash
ipython profile create yourprofilename

This has created a folder named profile_yourprofilename in your IPython folder. If you don’t know where this folder is, just type in a terminal

ipython locate

You have now a new profile ! In the future, to launch IPython with this new profile, use the command

ipython notebook --profile yourprofilename

To later rename a profile, go to your .ipython folder, rename the folder “profile_name” into “profile_newname”, then run

ipython profile create newname

For the moment your profile is identical to the default profile. Now go in your profile folder, and create a subfolder profile_yourprofile/static/css/. In this /css folder we will put CSS files and pictures to tune the appearance of the notebook. Start by creating a file custom.css, we will see how to fill it in the next sections. I just learned enough CSS to make the themes, so I am not really in control here. Comments and suggestions are welcome🙂

Theme 1: A Clean Notebook

My favorite🙂 . Here is the content of my custom.css . I just hope that the structure and class names of IPynb html pages will not change too much in the future

/* file custom.css for the theme CleanNotebook */
div#notebook { /* centers the page content */
padding-left:20%;
padding-right:20%;
}

div.cell { /* Tunes the space between cells */
margin-top:1em;
margin-bottom:1em;
}

div.text_cell_render h1 { /* Main titles bigger, centered */
font-size: 2.2em;
line-height:1.4em;
text-align:center;
}

div.text_cell_render h2 { /*  Parts names nearer from text */
margin-bottom: -0.4em;
}



div.text_cell_render { /* Customize text cells */
font-family: 'Times New Roman';
font-size:1.5em;
line-height:1.4em;
padding-left:3em;
padding-right:3em;
}

#notebook li { /* More space between bullet points */
margin-top:0.8em;
}

div.cell.code_cell {  /* Areat containing both code and output */
background-color:#F1F0FF; /* light blue */
border-radius: 10px; /* rounded borders = friendlier */
padding: 1em;
}

Bonus: Motivational Penguin

I am really fond of the motivational penguin from chibird.com (thanks Jacqueline !!!).
motivational_penguin

If you want the motivational penguin to cheer you up on the side of your screen while you are coding, put the GIF file in the /css folder and add these few lines to the code above:

#notebook_panel {
background : url('motivational_penguin.gif') no-repeat left center ;
}

Theme 2: Zen of IPython

For this theme you will need to place each of these pictures in your CSS folder (I made all of them from Public Domain pictures from Wikimedia Commons and I put all of them in the Public Domain, where they belong).

The css code is surely far from optimized.

div#notebook {
padding-left:20%;
padding-right:20%;
}

div#ipython-main-app {
background: url('Shotei.jpeg') no-repeat top left;
}

body {
background : url('paper.jpeg') ;
}

div.cell {
margin-top:20px;
margin-bottom:20px;
}


h1 {
text-align:center;
}


div#notebook li {
margin-bottom:20px;
}

p,li,h1 {
line-height:120%;
}

p,ul {
padding-left:10%;
padding-right:10%;
}

.text_cell_render{
font-family: 'Palatino','Arial';
font-size:1.8em;
}

.text_cell h1 {
font-size: 2.2em;
text-align:center;
line-height:120%;
}

.text_cell h2 {
font-size: 1.8em;
margin-bottom:0;
}


div.input_area {
background: url('mountain_paper.png') ;
border-radius: 15px;
overflow:auto;
}


div.CodeMirror {
background: url('mountain.jpeg') repeat-y top right;
background-size: 20%;
padding-right:20%;
font-size: 1.3em;
}

div.CodeMirror-scroll {
padding-bottom:1.01em;
}

.output_wrapper{
background-image: url('sky.jpeg');
border-radius: 15px;
}

div.output {
border-radius: 15px;
background: url('cherryFlower.jpeg') repeat-y right;
background-size: 25%;
}

div.prompt.input_prompt{
display:none;
}

div.code_cell {
background: #e6ddce;
-moz-border-radius: 15px;
border-radius: 15px;
padding: 15px;
font-size:14px;}

Theme 3: Gameboy

For this one you will need the buttons:
buttons

.prompt {display:none;}

div.cell.code_cell {
max-width: 800px;
background-color:white;
background:url('buttons.jpeg') no-repeat center bottom;
background-size:80%;
border-color: black;
border-radius: 30px;
padding-bottom:5px;
padding: 1em 1em 200px 1em;
}

div.output_wrapper {
margin-left:0px;
padding-left:0px;
padding-top:0px;
}

div.input_area, .output_stdout  {
margin-left:110px;
border-style:solid;
max-width:600px;
background-color: #626e02;
border-color: #656e7f;
border-width: 20px;

}

div.output_stdout {
border-top: none;
margin-top:-26px;
padding-top:26px;
}
div.CodeMirror-scroll {
padding-bottom:20px
}

div.CodeMirror, div.output_stdout pre {
padding: 0.5em 0.5em 0.5em 0.8em;
font-size: 1.6em;
font-family:"G.B.BOOT";
}

.cm-s-ipython { color: black; }
.cm-s-ipython span.cm-keyword {color: black;}
.cm-s-ipython span.cm-number {color: black;}
.cm-s-ipython span.cm-operator {color:black;}
.cm-s-ipython span.cm-meta {color: black;}
.cm-s-ipython span.cm-comment {color: black;}
.cm-s-ipython span.cm-string {color: black;}
.cm-s-ipython span.cm-error {color: black;}
.cm-s-ipython span.cm-builtin {color: black;}
.cm-s-ipython span.cm-variable {color: black;}

Delay differential equations in Python

I wrote a very simple and user-friendly method, that I called ddeint, to solve delay differential equations (DDEs) in Python, using the ODE solving capabilities of the Python package Scipy. As usual the code is available at the end of the post🙂.

Example 1 : Sine

I start with an example whose exact solution is known so that I can check that the algorithm works as expected. We consider the following equation:

y(t) = sin(t) \ \ \ for\ \ \ t < 0

y'(t) = y(t-3\pi/2) \ \ \ for\ \ \ t \geq 0

The trick here is that sin(t-3\pi/2) = cos(t)=sin'(t) so the exact solution of this equation is actually the sine function.

from pylab import *

model = lambda Y,t : Y(t - 3*pi/2) # Model
tt = linspace(0,50,10000) # Time start, end, and number of points/steps
g=sin # Expression of Y(t) before the integration interval 
yy = ddeint(model,g,tt) # Solving

# PLOTTING
fig,ax=subplots(1)
ax.plot(tt,yy,c='r',label="$y(t)$")
ax.plot(tt,sin(tt),c='b',label="$sin(t)$")
ax.set_ylim(ymax=2) # make room for the legend
ax.legend()
show()

dde_sin

The resulting plot compares our solution (red) with the exact solution (blue). See how our result eventually detaches itself from the actual solution as a consequence of many successive approximations ? As DDEs tend to create chaotic behaviors, you can expect the error to explode very fast. As I am no DDE expert, I would recommend checking for convergence in all cases, i.e. increasing the time resolution and see how it affects the result. Keep in mind that the past values of Y(t) are computed by interpolating the values of Y found at the previous integration points, so the more points you ask for, the more precise your result.

Example 2 : Delayed negative feedback

You can select the parameters of your model at integration time, like in Scipy’s ODE and odeint. As an example, imagine a product with degradation rate r, and whose production rate is negatively linked to the quantity of this same product at the time (t-d):

y(t) = 0 \ \ \ for\ \ \ t < 0

y'(t) = \dfrac{1}{1+(\dfrac{y(t-d)}{K})^2} -ry(t) \ \ \ for\ \ \ t \geq 0

We have three parameters that we can choose freely. For K = 0.1, d = 5, r = 1, we obtain oscillations !

from pylab import *

# MODEL, WITH UNKNOWN PARAMETERS
model = lambda Y,t,k,d,r :  1/(1+(Y(t-d)/k)**2) - r*Y(t)

# HISTORY
g = lambda t:0 

# SOLVING
tt = linspace(0,50,10000)
yy = ddeint(model,g,tt,fargs=( 0.1 , 5 , 1 )) # K = 0.1, d = 5, r = 1

# PLOTTING
fig,ax=subplots(1)
ax.plot(tt,yy,lw=2)
show()

dde_negativefeedback

Example 3 : Lotka-Volterra system with delay

The variable Y can be a vector, which means that you can solve DDE systems of several variables. Here is a version of the famous Lotka-Volterra two-variables system, where we introduce some delay d. For d=0 we find the solution of a classical Lotka-Volterra system, and for d non-nul, the system undergoes an important amplification:

\big(x(t), y(t)\big) = (1,2) \ \ \ for\ \ t < 0, \ \ else

x'(t) = 0.5x(t)\big(1-y(t-d)\big)\\   y'(t) = -0.5y(t)\big(1-x(t-d)\big)

from pylab import *

def model(Y,t,d):
    x,y = Y(t)
    xd,yd = Y(t-d)
    return array([0.5*x*(1-yd), -0.5*y*(1-xd)])

g = lambda t : array([1,2])
tt = linspace(2,30,20000)
fig,ax=subplots(1)

for d in [0, 0.2]:
    yy = ddeint(model,g,tt,fargs=(d,))
    # WE PLOT X AGAINST Y
    ax.plot(yy[:,0],yy[:,1],lw=2,label='delay = %.01f'%d)

ax.legend()
show()

dde_lotka

Example 4 : A DDE with varying delay

This time the delay depends on the value of Y(t) !

y(t) = 1,\ \ \ t \leq 0

y'(t) = - y\big(t-3\cos(y(t))^2 \big),\ \ \ t > 0

from pylab import *
model = lambda Y,t:  -Y(t-3*cos(Y(t))**2)
tt = linspace(0,30,2000)
yy = ddeint(model, lambda t:1, tt)
fig,ax=subplots(1)
ax.plot(tt,yy,lw=2)
show()

dde_Ydependent

Code

Explanations

The code is written on top of Scipy’s ‘ode’ (ordinary differential equation) class, which accepts differential equations under the form

model(Y,t) = ” expression of Y'(t) ”

where Y, and the output Y', must be Numpy arrays (i.e. vectors).

For our needs, we need the input Y to be a function of time, more precisely a function that can compute Y(t) at any past or present t using the values of Y already computed. We also need Y(t) to return the value of some function g(t) if t is inferior to some time tc that marks the start of the integration.

To this end, I first implemented a class (ddeVar) of variables/functions which can be called at any time t: for $t<latex t_c$, it will return the value of g(t), and for t>tc, it will look for two already computed values Y_a and Y_b at times t_a<t<t_b, from which it will deduce Y(t) using a linear interpolation. Scipy offers many other kinds of interpolation, but these will be slower and won't support vectors for Y.

Such variables need to be updated every time a new value of Y(t) is computed, so I created a class 'dde' that inherits from Scipy's 'ode' class but overwrites its integration method so that our special function Y is updated after each integration step. Since 'ode' would feed the model with a vector Y (a Numpy array to be precise), which we don't want, we give to the integrator an interface function that takes a Numpy array Y as an argument, but immediately dumps it and calls the model with our special ddeVar variable Y (I hope that was clear🙂 ).

Ok, here you are for the code

You will find the code and all the examples as an IPython notebook HERE (if you are a scientific pythonist and you don’t know about the IPython notebook, you are really missing something !). Just change the extension to .ipynb to be able to open it. In case you just asked for the code:

# REQUIRES PACKAGES Numpy AND Scipy INSTALLED
import numpy as np
import scipy.integrate
import scipy.interpolate
 
class ddeVar:
    """ special function-like variables for the integration of DDEs """
    
    
    def __init__(self,g,tc=0):
        """ g(t) = expression of Y(t) for t<tc """
        
        self.g = g
        self.tc= tc
        # We must fill the interpolator with 2 points minimum
        self.itpr = scipy.interpolate.interp1d(
            np.array([tc-1,tc]), # X
            np.array([self.g(tc),self.g(tc)]).T, # Y
            kind='linear', bounds_error=False,
            fill_value = self.g(tc))
            
            
    def update(self,t,Y):
        """ Add one new (ti,yi) to the interpolator """
        
        self.itpr.x = np.hstack([self.itpr.x, [t]])
        Y2 = Y if (Y.size==1) else np.array([Y]).T
        self.itpr.y = np.hstack([self.itpr.y, Y2])
        self.itpr.fill_value = Y
        
        
    def __call__(self,t=0):
        """ Y(t) will return the instance's value at time t """
        
        return (self.g(t) if (t<=self.tc) else self.itpr(t))
 


class dde(scipy.integrate.ode):
    """ Overwrites a few functions of scipy.integrate.ode"""
    
    
    def __init__(self,f,jac=None):
		
        def f2(t,y,args):
            return f(self.Y,t,*args)
        scipy.integrate.ode.__init__(self,f2,jac)
        self.set_f_params(None)
        
 
    def integrate(self, t, step=0, relax=0):
		
        scipy.integrate.ode.integrate(self,t,step,relax)
        self.Y.update(self.t,self.y)
        return self.y
        
 
    def set_initial_value(self,Y):
		
        self.Y = Y #!!! Y will be modified during integration
        scipy.integrate.ode.set_initial_value(self, Y(Y.tc), Y.tc)
 
 
 
def ddeint(func,g,tt,fargs=None):
    """ similar to scipy.integrate.odeint. Solves the DDE system
        defined by func at the times tt with 'history function' g
        and potential additional arguments for the model, fargs
    """
    
    dde_ = dde(func)
    dde_.set_initial_value(ddeVar(g,tt[0]))
    dde_.set_f_params(fargs if fargs else [])
    return np.array([g(tt[0])]+[dde_.integrate(dde_.t + dt)
                                 for dt in np.diff(tt)])

Other implementations

If you need a faster or more reliable implementation, have a look at the packages pyDDE and pydelay, which seem both very serious but are less friendly in their syntax.

Typing keyboard + python = Musical instrument !

In this post I will show you how to do this :

I am not the first to do that, but most people who play their computer on Youtube use either very expensive programs, or programs that won’t run on your computer, or not-so-efficient programs with not-that-much possibilities of extension, or cheap programs with a big lag between pressing the key and actually hearing the note.

So here is a very small Python script which will run fine even on a basic netbook. If you are not familiar with Python, you should take online courses , it is really worth it🙂 !
If you are faminiliar with Python, then you are welcome to improve the code on its Github page.

Transforming your keyboard into a mixtable

My original idea was to transform my computer keyboard into a mixtable, to make something like in this very awesome video:

So my first move was to make a program that would take a configuration file my_configuration.cf containing this:


q, dog.wav
w, cat.wav
e, tarzan.wav
r, scream.mp3

And then if you hit q you would hear a dog from the dog.wav soundfile, if you hit w you’d hear a cat, etc…
This is pretty easy to do with Python’s pygame package. Here is my code (inspired by similar stuff from the package Mingus):


import pygame as pg
import csv
import time


SAMPLE_WIDTH = 16
FPS = 44100
N_CHANNELS = 2
BUFFER = 2**9
    
def minimix(config_file,mode = 'instrument'):
    """
    Opens an interface that lets you press the keys of your keyboard
    to plays the related soundfiles as sepcified in the provided
    configuration file.
    
    Args:
       config_file (str):  A file associating keyboard keys with
                           file names. Each line must be of the form
                           key , filename.
                       
       mode (str) :
            instrument -- causes autorepeat of samples as long
                            as the key is pressed.
            sustain -- causes autorepeat of the samples until its
                       key is pressed again.
            player -- striking the key causes the sound to play once. 
    
    Returns:
       a list of the (time,events).
    """
    
    repeat = 0 if (mode is 'player') else (-1)
    
    pg.mixer.pre_init(FPS,-SAMPLE_WIDTH,N_CHANNELS,BUFFER)
    pg.init()
    screen = pg.display.set_mode((640,480))
    
    
    
    
    ##### READ CONF FILE
    
    key2sound = {}
    key2file = {}
    config = csv.reader(open(config_file, 'rb'), delimiter=',')
    
    for key, soundfile in config:
        
        key,soundfile = key.strip(' '),soundfile.strip(' ')
        
        if key is not '#':
            
            key2file[key] = soundfile
            key2sound[key] = pg.mixer.Sound(soundfile)
    
    events_list = []
    currently_playing = {k : False for k in key2sound.iterkeys()}
    
    
    ##### MAIN LOOP
    
    while True:
        
        event =  pg.event.wait()
      
        if event.type in (pg.KEYDOWN,pg.KEYUP):
            key = pg.key.name(event.key)
              
            if key in key2sound:
                
                if event.type == pg.KEYDOWN:
                    
                    if (mode == 'sustain') and currently_playing[key]:
                        
                        key2sound[key].fadeout(20)
                        currently_playing[key] = False
                        
                    else:
                        
                        key2sound[key].play(repeat) 
                        currently_playing[key] = True
                    
                    events_list.append((time.time(),key2file[key]))
                    
                elif event.type == pg.KEYUP and (mode == 'instrument'):
                    
                    key2sound[key].stop()
                    currently_playing[key] = False
                    
                    events_list.append((time.time(),key2file[key]))
            
            elif event.key == pg.K_ESCAPE:
                
                break
    
    pg.quit()
    
    return events_list

Transforming your keyboard into a musical instrument

If instead of using various noises like cat and dog you use different notes from the same instrument, then you turned your computer into some kind of piano. The problem is that a set of soundfiles with all the notes of an instrument is difficult to find on the internet, so I wrote a script that makes as many notes as you want from just one sound by shifting its pitch up or down. It uses the audio processing program Soundstretch that you will need to install first :

import os

def shift_wav(wavfile,output,shifts,verbose=False):
    """
    Makes new sounds by shifting the pitch of a sound.
    Requires soundstrech installed.
    
    Args:
        wavfile : Name of the file containing the original sound
        output: name to use as a prefix for the output files and for
                the output folder name
        shifts (list of int): specifies of how many half-tones the pitch
                shifts should be. For instance [-2,-1,1,2] will produce
                4 files containing the sound 2 half-tones lower, one
                halftone lower, one halftone higher and two halftones
                higher.
    """
    
    folder = os.path.dirname(output)
        
    if not os.path.exists(folder):
        
        os.makedirs(folder)
        
    for i,s in enumerate(shifts):
        
        outputfile = '%s%02d.wav'%(output,i)
        
        command = 'soundstretch %s %s -pitch=%d'%(wavfile,outputfile,s)
        if verbose:
            print command
        os.system(command)

Going further

There is so much one could do to improve the program.

On the musical side, for instance, finding configurations of the keyboard that are particularly ergonomic. In the video above I used this configuration:

azerty

I called it typewriter because it enables you to play very fast things while moving your hands a minimum (the video is a bad example :P) . But maybe there is better to find !

Also, one could start listing every cool piece of music that can be played on a typing keyboard. I use 46 keys ( almost 4 octaves !), that makes a lot of possibilities !

On the programming side, there is a lot of little things I can think of, like automatizing scale changes, introducing nuances, designing nice interfaces (why not a guitar-hero-like game where you would actually be playing music on a playback ?), writing a script that would take some sheet music (in a nice format, like ABC, MIDI, lylipond) and return the list of the keys you should strike to play it.

I actually wrote a lot more code, for instance to make it easier to write configuration files, for sound processing, etc., but since it is not strictly necessary I am not reporting it here (I’ll certainly put a working version on GitHub, or such, some day).

Philosophy of the musical keyboard

Your typing keyboard is a real instrument. Of course it is not its primary use, but our voice’s primary purpose was not to sing, either. Now do the math : how many people out there have a piano at home ? And how many have a computer ? That gives you an idea of how many people would like to play the piano, cannot, but could play their computers instead.

So promoting computer-keyboardism is ultimately about bringing music to the masses. It is about providing everyone with an instrument that you can practice at home, in the train, at work, and that will be familiar to anyone everywhere in the world.

There is more : how many of you, pianist readers, have started the piano for seduction purposes ? (yeah, sure, me neither…) But public places with a piano on which you could show off your mad skills are getting pretty rare, aren’t they ? Especially since most bars have traded their good old piano for a TV. But think about all these places with a computer at hand ! Yep, time for you to develop a talent that will be useful in real life !

So practice, get good, be one of the first composers for tomorrow’s instrument, impress your friends and spread the good news ! If you are still reading me after so much gibberish , then do not hesitate : you just proved how little you value your free time, you are the right person for the task !

An analytical solution to a Jewish problem

Maybe you have heard about the article of Khovanova and Radul, Jewish Problems, in which they collect tricky mathemical problems that were alledgedly designed to prevent Jews (to which they were specifically given) from passing the oral entrance exams of Moscow State University . With such a story and such a naïvely equivocal title, you can count on all the comments-section jewishologists of the world to promote your article🙂

What I found interesting is this article is that, of the 21 problems given, I could answer most of the analytical ones, while I was completely unequipped to solve the geometrical problems. Has geometry always been that hard ? My impression is that, at least here in France, its teaching is slowly being replaced by linear algebra . Maybe it is simply less needed today. Take civil engineering : a beam used to be some kind of parralelepipede, now it is a grid of finite elements !

As an illustration, let us have a look at the fifth problem: solve the equation

\sin^7(x) + \dfrac{1}{\sin^3x} = \cos^7(x) + \dfrac{1}{\cos^3x}.

What kind of math problem do you see here ? (let’s just hope it has nothing to do with geometry !) The authors of the article give a solution based on some transformation of the equations and a few trigonometricks (at some point they use the variable t=sin(x)cos(x) and retransform it into sin(2x)/2). What I first saw when I read the question (and I bet that’s what you saw too if you have my kind of training) was an analytical problem involving the function

f(t) = t^7 + \dfrac{1}{t^3}

for which it is asked to find all the (u,v) which are the sine and cosine of a same angle and verify f(u)=f(v) ! An obvious situation where this will work is when u=v, which leads us to the solutions t = \pi/4 and t = -3\pi/4, which are the only angles whose sine is equal to the cosine. Now, are there any other solutions ? In the end what is asked is whether the function f is bijective, or not, or just enough ! Let us compute its derivative:

f'(t) = 7t^6 - 3 \dfrac{1}{t^4} = \dfrac{7t^{10}-3}{t^4}

This is only positive outside the interval [-\sqrt[10]{\frac{7}{3}} , +\sqrt[10]{\frac{7}{3}}] . We can now sketch the function:

Since u and v are sine and cosine, they will belong to the interval [-1,1]. In red I have represented two non-bijectivity zones of the function on this interval: if u and v are two different numbers of [-1,1] verifying f(u)=f(v), then the two of them must belong to one of these red zones. Now, notice that
f(\sqrt{\frac{1}{2}}) = (\sqrt{\frac{1}{2}})^7 + \sqrt{2}^3 > \sqrt{2}^3 = 2\sqrt{2} > 2
This shows that \sqrt{\frac{1}{2}} is placed left to the red zone in the right (see figure). A consequence is that, if u and v belong to the right red zone, then

u^2 + v^2 > \sqrt{\frac{1}{2}}^2 + \sqrt{\frac{1}{2}}^2 = 1

so u and v cannot be the sine and the cosine of the same angle (or their squares would sum up to 1). By imparity of the function, the same can be said about the left red zone. As a conclusion, it is not possible that f(u)=f(v) if u and v are the sine and cosine of a same angle and have different values. So t = \pi/4 and t = -3\pi/4 are the only two solutions to the given equation. Did anyone see another solution ?

Counting bacteria : confidence intervals from one measurement

Here are two biologically-inspired statistics problem, with very simple solutions. They show that, when you have informations about the way the observations are collected, one observation may be enough to give a confidence interval.

Question 1

From a bacterial solution I sampled a small volume v in which I found no bacteria. Give an upper bound (with confidence 95%) for the concentration of bacteria in the solution.

Question 2

From a bacterial solution I sampled a small volume v in which I found n bacteria (n > 20). Give a 95% confidence interval for the concentration of bacteria in the solution.

Solutions

Short answers : for the first question c < (3/v) and for the second c \in [\frac{ n \pm 2\sqrt{n}}{v}] .

Yep, it’s that simple ! See the solutions.

Physics of Harrison’s H1 clock

Today’s post is about a subject I presented a few years back at the entrance exams of the french engineering schools. If you like clocks, boats, adventures and unnecessary complications, this is for you🙂

Historical notice

During the Renaissance the English realised how good they were on seas, and how much potential there was for some ass-kicking with all these loaded European ships making round trips to America. A very crucial problem in sailing is to know where you are. For that you need

  1. Your latitude : how far north are you from the equator ?
  2. Your longitude : how far west are you from London ?

Getting the latitude is easy:

  • During the day, look at the sun (especially at noon). The higher it is, the nearer you are from the equator
  • At night, look at a pole star : the higher it is, the farther away you are from the equator.

But no such tricks to estimate the longitude existed in the 1700’s. Astronomers had thought of something but it was all very complicated and stuff. So the English did what you do when you are out of ideas: they made a law (the Longitude Act) proposing a huge prize (that they had no intention of giving up too easily) to any subject of his Majesty who would come up first with a method.
That’s were John Harrison comes in, with a beautifully simple idea:

  1. Take a clock indicating London’s time onboard.
  2. During your travel, look at the sun to get the time of the place where you are, and compare it to the clock’s time.
  3. If for instance the sun indicates 12h00, and the clock indicates that it is 10h00 in London, then you know that you are \frac{12-10}{24}*360 = 30 degrees west from London.

The only problem is that, if you take your good old pendulum clock on a boat, it will instantly go mad due to all that pitching, rolling and yawing ! So you need to design a more sophisticated device. But John Harrison had a slight advantage here: he was a professional clockmaker !

John Harrison (1693-1776) Credit: Wikipedia

Shake it baby ! Harrison’s H1 clock

In 1730, Harrison imagines his first marine clock, that he calls H1:

Harrison’s H1. Credit : National Maritime Museum, (Greenwich)

Here is an idealized scheme of the beast’s arms :

I believe that the contact between the two arms is there to make sure that they will always move in symetrical ways (otherwise the movement could get chaotic). A very nice feature of this design is that any force that is applied equally to each of the weights is cancelled out and has no effect on the oscillations !

This makes the clock’s movement independent from gravity (meaning that the clock doesn’t need to be vertical to oscillate properly) and from any translational movement of the boat !
However (and that’s where my study started) I noticed that the centrifuge force induced by the rolling of the ship is not compensated, as the weights are not equally distant from the rotation center:

And my (very fundamental and of much historical importance) question was : how much does the rolling affect the clock ?

A small study

Before getting to the real question, let us study how this clock works when it is not perturbated.
Suppose that the left arm, of radius R makes an angle \theta with its equilibrium position (that we suppose to be vertical):

The other arm is in a symetrical position, so the upper spring is streched by a length 2R\sin(\theta) and thus pulls the upper weight of the left arm with a force equal to 2kR\sin(\theta) (where k is the spring’s rigidity constant). The lower spring is compressed by a length 2R\sin(\theta) and the lower weight undergoes a push of 2kR\sin(\theta). These two forces apply negative torques on the arm, both with a lever equal to R\cos(\theta). Finally, neglecting the mass of the arm’s sticks, the arm’s angular momentum is 2R^2M\dot{\theta} where M is the mass of one weight. Once you have all that, the angular momentum theorem gives

2R^2M\ddot{\theta} = -2Rk\sin(\theta)R\cos(\theta)

which after simplifications gives

\ddot{\theta} + \dfrac{k}{M}\sin(2\theta) = 0

Note that (partly because my model is very simple) the oscillations do not depend on the radius of the arms. If the amplitude of the oscillations if really small, we have \sin(2\theta) \simeq 2\theta and the equation becomes

\ddot{\theta} + 2\dfrac{k}{M}\theta = 0

where we recognize an harmonic oscillator of period T = \pi \sqrt{\dfrac{2M}{k}} . If the oscillations are bigger, one can still simulate the dynamics (here with Python):


from pylab import *
from scipy.integrate import odeint

def H1diff(Y,t,k,M):
    """
    Equations describing the dynamics of the model
    Y is [A,dA/dt], t is the current time, k is the spring's constant,
    M is the mass of the weights.
    Returns [dA/dt, ddA/ddt] according to
    dd(A)/ddt + k/M sin(2A) = 0
    """

    A,dA = Y
    return array([dA, -k / M * sin(2*A)])

k=1.0 ; M=1.0 ; initialAngle = pi/6 # parameters
Y0= [pi/6, 0.0] # initial angle pi/6, initial speed 0
ts=linspace(0,10, 100) # 100 times points between 0 and 10
angles = odeint(H1diff,[initialAngle,0],ts,args=(k,M))[:,0] # solving

# plotting
fig,ax = subplots(1)
ax.plot(ts,angles) # plotting
ax.set_xlabel('time',fontsize = 24)
ax.set_ylabel('angle',fontsize = 24)
show()

Or better🙂


def makeAnimation(k=1.0,M=1.0,nframes=40,R=1.0,**kwargs):
    """ simulates and creates a GIF of the H1 """

    # default for tmax is one period.
    tmax = kwargs.pop('tmax',2*pi*sqrt(M/k))

    # delay is computed so that the animation actually lasts tmax.
    delay = kwargs.pop('delay',int(100.0*tmax/nframes))

    O1 = array([-1.0,0]) # center of the left pendulum
    O2 = array([1.0,0])  # center of the right pendulum
    R = 1.0

    def drawH1(angle,ax):
        """ draws an Harrison H1 clock given the angle of the arms """

        # compute the coordinates of the four weights:

        s,c = sin(angle), cos(angle)
        LU = O1 + R*array([-s,c])
        LD = O1 + R*array([s,-c])
        RU = O2 + R*array([s,c])
        RD = O2 + R*array([-s,-c])

        # draw the springs

        for w1,w2 in (LU,RU),(LD,RD):
            ax.add_line(Line2D(*zip(w1,w2), linewidth=2,color='b'))

        # draw the sticks

        for w1,w2 in (LU,LD),(RU,RD):
            ax.add_line(Line2D(*zip(w1,w2), linewidth=4,color='k'))

        # draw the joints

        for point in O1,O2:
            ax.add_artist(Circle(point,radius = 0.1, color='grey'))

        # draw the weights

        for point in LU,LD,RU,RD:
            ax.add_artist(Circle(point,radius = 0.2,color='k'))

    Y0 = [pi/6,0]
    ts=linspace(0,tmax,nframes)
    angles = odeint(H1diff,Y0,ts,args=(k,M))[:,0].flatten()

    fig,ax = subplots(1)
    fig.set_size_inches(2.5,2.5)
    ax.set_xlim(-1.8,1.8)
    ax.set_ylim(-1.3,1.3)

    for i,angle in enumerate(angles):
        ax.clear()
        axis('off')
        drawH1(angle,ax)
        fig.savefig('hhh%.3d.jpeg'%i)

    # make the gif

    os.system('convert -delay %f -loop -1 hhh*.jpeg harrison.gif'%delay)

    # remove the temporary files
    for myFile in os.listdir('./'):
        if myFile.startswith('hhh'):
            os.remove(myFile)

H1 simulation

Any role for the Rolling ?

Back to our previous question: how much does the rolling affect the clock’s movement ? Let us just put a few things straight first:

  • It’s impossible to come even close to a solution that I could defend more than one minute.
  • Anyway, who cares ! It’s not like if it was an important problem.

This being said, here is how I went with the modelling. For the clock:

  1. I emailed the Greenwich Museum, where the original H1 is, and they were nice enough to give me a few measures of the clock (length of the arms, radius and material of the weights)
  2. I estimated the maximal angle and checked the period using a video found on the internet.
  3. I then determined the constant k of the springs by simulation and optimisation, as the optimal constant to get a period of one second.

Here is how you do the last step with Python:

from scipy.optimize import fmin

def find_k(T,M,initialAngle):
    """ finds the optimal spring rigidity constant k to obtain an H1
        clock with a period of T given the mass M of the weights
        and the initial angle of the clock """

    def f(k):
        ts=[0,T]
        res = odeint(H1diff,[initialAngle,0],ts,args=(k,M))
        finalAngle = res[-1,0]
        return (finalAngle-initialAngle)**2

    k_init = 2*(pi**2)*M / T**2  # expected for very small oscillations
    return fmin( f ,k_init)

To model the boat:

    1. I found a plan on the internet of the very famous HMS Bounty, which is approximately the kind of ship used in Harrison’s times
  1. Since the plan was accurate enough to make a mesh out of it, I used pens and rulers to get the coordinates of the points of the slices of the ship on the plan. then I used some sort of finite elements method (I had not heard of this at the time) to simulate the oscillations of the boat on still water. This gave me a profile of the ship’s rolling over time.
  2. A cool thing : I found on the internet that some people in Australia had made an exact replicate of the Bounty. So I emailed the captain to see if he could measure how long the ship takes to do one roll. He found exactly the same period, more or less one second !

Putting everything together I was able to see how much the rolling impacts the period of the clock over a few minutes, then how much difference it makes in the estimated time after a few days, and how much this represents in terms of distance… Results that I don’t really remember, but once again, who cares ? Isn’t the most important having fun, squeezing equations and poking people from all around the world ?

Epilogue (with spoilers)

The H1 wasn’t precise enough to compute the longitude as precisely as demanded by the Longitude Act. Further improvements kept Harrison busy for his whole life. He produced three clocks, H1, H2, and H3, until he realised that he had been thinking too big, and that the solution might come from a bleeding edge technology of the time: the chronometer ! With its small size and it’s reliable spiral spring, this ancester of the pocket watch was accurate enough to finally get Harrison his prize, after more than 40 years of research !