The Model

This project is to write a program to work out the phonons (normal mode vibrations) of a Fibonacci quasicrystal with two different sorts of atom. The ``adults'' and ``children'' are the masses, $M$ and $m$, of the two kinds of atom. Since we are only interested in small vibrations, we may represent the forces between the atoms as simple springs (spring constant $K$) and we will assume that all the springs are identical ($K$ is independent of the types of the atoms at the ends of the spring). The equation of motion of such a system may be written as

$$m_n {\d^2 x_n\over\d t^2} = -m_n\omega^2 x_n = K\left(x_{n+1} - 2 x_n + x_{n-1}\right) \quad,$$ (3.43)

where $x_n$ and $m_n$ are the displacement and mass of the $n$th atom. If the chain contains $N$ atoms, you will then have a set of $N$ coupled second order ordinary differential equations.

The choice of boundary conditions for the atoms at either end of the chain is up to you and should not make much difference when the chain is long enough: you could fix the atoms at either end to immovable walls, $x_0 = x_{N+1} = 0$, or you could leave them free by removing the springs at the 2 ends, $K_{0,1} = K_{N,N+1} = 0$. (Periodic boundary conditions -- when the chain of atoms is looped around and joined up in a ring -- are convenient for analytic work but not so good for numerical work in this case. Why?)

The equations (3.43) are $N$ linear algebraic equations which may be cast as a tridiagonal matrix eigenproblem. The eigenvalues of this problem are $\omega^{2}$ and so give the vibrational frequencies, and the eigenvectors give the corresponding normal mode coordinates, $x_{n}, \quad n=1,N$.

Note, however, that the presence of the masses $m_n$ in (3.43) means that the problem is in the generalised form $\bss{A}\bi{x} = \omega^2 \bss{B}\bi{x}$. As discussed in the notes and in problem 6 this can be transformed into the normal eigenvalue problem by making the substitution $x_n = m_n^{-{\mathchoice{{\textstyle{\frac12}}}{{\textstyle{\frac12}}}{{\scriptstyle{1/2}}}{{\scriptscriptstyle{1/2}}}}} y_n$ and multiplying the $n$th equation by $m_n^{-{\mathchoice{{\textstyle{\frac12}}}{{\textstyle{\frac12}}}{{\scriptstyle{1/2}}}{{\scriptscriptstyle{1/2}}}}}$. I suggest you solve the eigenproblem by using a NAG or LaPack routine for the eigenvalues only. There should be no problem in dealing with chains of several hundred atoms or more.