Spring Time

Picture a spring. You know that the harder you squash or stretch it the harder it tries to get back to its "natural length". So you know that the force must be related to the difference between the natural length and the length it's been forced it to. \[F = f(dx)\] where $dx$ is the difference between the natural length and the stretched or squashed length. The equation above summarizes this in mathematical form, the force is given by some function of $dx$. $17^{\text{th}}$ century weirdo Robert Hooke gave the first guess for what this function was in the form of a Latin anagram ceiiinosssttuv which is solved to give Ut tensio, sic vis, "As the extension, so the force" or "Force is proportional to the extension." This gives us \[F = k dx \] where $k$ is some constant that tells us how stiff the spring is. Of course as advanced $21^{\text{st}}$ century scientists we can put in a more complicated function, \[F = k dx + \alpha dx^2 + \beta dx^3 \] $\alpha$ and $\beta$ are just some numbers. Any force law, other than Hooke's one is called "non-linear".

Now you could take a bunch of springs and join them together, making something like a mattress, pull on it and let it oscillate. For springs obeying Hooke's law it will oscillate in a pretty boring way. For non-linear springs the motion will be crazier. Joining a bunch of non-linear springs together is called the Fermi-Pasta-Ulam problem. Fermi and Ulam are famous guys. Pasta was a American computer scientist with an unfortunate name. Also involved in the first paper was a lady called Mary Tsingou, but for some reason she doesn't get to be included with the first three guys in the name of the problem. The idea was that if you had a non-linear force the system would oscillate in such a crazy way that it would go through every possible arrangement. When they actually simulated it however they found it didn't, it oscillated in a surprisingly regular way. The conclusion in fancy physics jargon: non-linearity does not imply ergodicity.

From the force we can get a potential energy, using the non-linear force above gives us the following four basic shapes of the potential.

You can think of the system as a little ball rolling on the potential hill. The top left is Hooke's law, small oscillations about the mean value, imagine a little ball rolling back and forth inside the "well". If you go to the simulation below you can see what that looks like by setting $\alpha = \beta = 0$ and $k = 1$ (and clicking Submit) for example. Much more interesting is the top right potential. Imagine carefully placing a little ball on top of the hump in the middle. Unless you were extremely careful it would roll down to the left or right. This is in fact exactly the sort of thing that happens in the Higgs mechanism which explains the mass of the electron. In our model, if we start almost at the top of the hump (press the "random" button and try $k=-1$, $\beta = 80$, $\alpha = 0$) we find that the springs want to compress together making clumps. These clumps are a little like our Higgs bosons. The symmetry, where all the points were more or less equal, is broken and we rolled down the hill.

The third one is also interesting, again imagine a little ball rolling on the hill on the bottom left. If it starts near the middle it will roll happily inside the little well, but if it gets above the hill it will roll down forever getting faster and faster. You can see this happen, and break the simulation, by choosing say $k = 1$ and $\beta = -250$. After a little while some of the springs should, by chance, accumulate enough energy to get over the hump and fall to their doom (infinity). The bottom right potential is what it looks like if you turn $\beta = 0$ and turn on $\alpha$. You have one stable direction and one unstable one.

k:

α:

β:

Try clicking on the dots, you can move them around if you want a different starting position. You can set up some interesting looking oscillations by choosing a good starting point for all the points. If you really want to, you can break some springs. Easiest is to make a potential with positive parameters. The forces can become huge and the algorithm that evolves the positions forward in time will break as the numbers overflow the maximum allowed size. These masses will break off. By pulling really hard on a spring you might be able to get it to snap off also, depending on the parameters. If you do break some springs clicking "random" will regenerate them.

A Revolutionary Model

As another application of an Ising type system in 'social physics' I am going to look at the paper, Peer to Peer and Mass Communication Effect on Revolution Dynamics, if you want more details you can read it, I have nothing to do with the authors and my mistakes are my own. The essential idea is again to simplify people down to binary agents with two possible choices, join the revolution or stay with the status quo. The agents have a predetermined conservatism and they can be influenced by their immediate neighbours. In addition the revolutionary agents have a youtube channel, a 50 megawatt broadcasting tower or a guy handing out leaflets. They have some way to make their influence felt by all the other agents simultaneously.

As an Ising system conservatism is something like a local coupling constant, the neighbour interaction is the usual ising interaction and the mass communication is like a global field, but created by the spins themselves. These analogies are far from exact however and as usual the best way to find out what the model does is to simulate it. We start with an initial state where most people are in the status quo and there are a few randomly distributed revolutionaries. Like the Schelling model but unlike the Ising model we don't allow people to change their minds and go back to the status quo. We could use the standard Metropolis accept reject but it wouldn't modify the conclusions much, as it is you can think of this as the Metropolis algorithm running at zero temperature.

a:

J:

c:

% rev:

Let $R = N_{rev}/N$ be the ratio of revolutionaries to the total population, initially $R = 0.1$. The parameter $a$ controls the power the mass communication has to influence people: given R let $H(R) = a R^b$ this is the 'media field'. I just follow the paper and set $b = 0.2$ so $a$ controls the strength of the mass media. $J$ controls the influence of the local environment, given an agent $i$ let $H(i) = J \sum_j S_j$ where the sum is over neighbours of $i$, in this case up down left and right. The final parameter $c$ controls the conservatism. Let $h_i$ be the conservatism of the agent $i$. Low $h_i$ means the agent doesn't much like the status quo while an agent with large $h_i$ is very attached to it. The $h_i$'s are randomly distributed at the start of the simulation (whenever you click the random button or move the $c$ slider) according to $P(h_i < h) = h^c$ for $h \in [0,1]$. For low c this gives a wide distribution of people with varied opinions about the current establishment. For large $c$, everyone is about equally conservative. Finally an agent joins the revolution if at any time, \[ h_i + H(i) < H\]

There are several regimes. Let's take $c = 2$ first. When nearest neighbours have little influence, $J \approx 0$ and the effects of mass media are also relatively small $a < 1$, we see an inhomogeneous distribution, much like a high temperature Ising model. There is no clustering of revolutionaries with revolutionaries and conservatives with conservatives. This is probably a dangerous situation in a real society, maybe this is a civil war phase. As the power of the media increases there is total revolution. For moderate $a$, say $0.3$ and moderate $J$, say $0.25$, we start to see clustering. Now that the influence of neighbours becomes as important as the influence of the media we start to see small clusters of conservatives swallowed up by the revolution, but larger clusters surviving. Perhaps we end up with two independent nations in this phase. Unsurprisingly if you increase the conservatism it requires stronger media influence and stronger neighbour coupling to get total revolutions. There also appears to be no phase where there are large numbers of revolutionaries and large clusters of conservatives coexisting. So no civil wars and no new nations, maybe not surprising considering a high $c$ nation is a nation where everyone is almost equally conservative.

People keep finding new applications for this old model which seems apt to describe situations where there is clustering, as here and in the Schelling model, or where there is a change of phase as here in the total revolution case and in the original Ising model at $T_c$. Where there is a binary choice; spin up or down, move or stay, revolt or don't, an Ising like model seems to often reproduce much of the interesting behaviour.

The Schelling Model

Lest you think I've forgotten about the Ising model, here is another application (sort of). Below are some maps of cities in America where each dot is a person and the colour represents their race. Taken from here. As you can see the distribution is less than homogeneous. I especially like the little red enclave in Detroit as well at the obvious "domain wall" along 8-Mile Road. New York and Chicago have more types of people. New York has many smaller clusters, possibly helped to form by geography.

New York

Chicago

Detroit

An economist named Thomas Schelling heard about these divisions and wanted to know if they were caused by people having strong preferences to be near their own kind, so he created a model to get to the essence of the problem. Like any good model we strip away the inessentials. People move for various reasons; for new jobs, better schools, to avoid meteor strikes, in Schelling's model people move based only on who their neighbours are. Our city is a large grid. At each point on the grid we have an agent of type A, an agent of type B or an empty space. When we examine the grid points we count the total number of neighbours and the number of neighbours of the same type as the central point. Note we are looking at the Moore neighbourhood of A.

In the plot above A has 7 neighbours 3 of which are the same type as him and 4 of which are different. We now pick a number called the intolerance which acts very much like the temperature in the Ising model. If the equation \[ \frac{ \text{Number different neighbours} }{\text{Number of neighbours}} > \text{intolerance} \] is true the agent moves to a vacant space (chosen at random), if not the agent stays put. For example with an intolerance of $30\%$ the agent will be okay with being in the minority as long as at least $30\%$ of its neighbours are the same type as it. Click the box below to see this in action, starting from a random distribution of agents and evolving.

intolerance:

occupancy:

Types:

% same:

The process we sought to model, racial segregation, has materialised. The original somewhat surprising finding was that for quite tolerant agents the final "diversity" measured by checking how many neighbours of the same type each agent has is relatively low (or the % same = 1 - diversity is relatively high). Much lower than what each individual agent will tolerate. Slightly biased individuals give rise to a very biased society. As usual there are some parameters to play with. Change the intolerance to see a phase transition to a noisy phase where everyone hates being near agents of different type so much that they move all the time. Decrease the occupancy to make the grid more sparsely populated, you can see that the segregated domains are separated by a "wall" of vacant sites. Finally play with the number of agent types to see that even an ethnically diverse city segregates, though in this case the place looks a bit more diverse since there are more, smaller clusters, kind of like the New York map. Interestingly, if the segregation phenomenon occurs then the agents are made extremely tolerant this does nothing to increase diversity, proving mathematically that society can never improve :)

Going back to look at the Ising model the differences are (a) conservation of agents: no one dies during this simulation - there are always the same number of A type and B types. The Ising model doesn't conserve number of up spins for example though it will fluctuate pretty close to its equilibrium value. (b) However this means there is no "conquest" of A agents by B agents. For example if you put the Ising model to a random configuration and rapidly lower the temperature (increase $\beta$) you will see a similar domain structure form. However, wait long enough and one spin will ultimately dominate. I can't decide if this makes a more or less realistic model. Also if you play with the Ising model in this way you sometimes get two solid blocks of up and down separated by a domain wall, Detroit style, which I haven't seen occur in the Schelling model.

Playing Tetris is my homework!

One summer when I was younger, twelve or so, I was sent to a summer school for the insufferably gifted. I really didn't want to go so I refused to pick a class, which was a terrible idea. In the end I was sent to learn computer programming, something I never would have chosen, since I liked video games. Being denied a summer of watching daytime TV and playing Street Fighter was an outrage and in protest I slept through most of the classes.

There was a final project that we had to do, based on what we learned during the course. People did very simple things like programming a computer to deal blackjack or do a substitution cypher. Things like that. But one kid, and I wonder what happened to him, made a Tetris program. This blew my mind. So, almost 15 years late, I am handing in my homework. I promise I didn't copy off anyone.

Click the left hand screen to begin playing. The left and right arrow keys move the blocks, down arrow slams them. The 's' and 'd' keys rotate the blocks. The 'p' key pauses and unpauses it. If you die click the left screen again to restart. Tetris is such a classic it can hardly be improved, however if you want to try something funny click on the box that says "Fixed" and try playing Tetorus.

The Doodle of Arthur

I think there's a lot in a name, probably because I happen to have an unusual one. Maths seems to yield topics and titles that for some reason always sound interesting to me. A random sample of some recent papers, about which I have no idea: "Non-Abelian Lie algebroids over jet spaces", "The Voevodsky motive of a rank one semiabelian variety", "On Avoiding Sufficiently Long Abelian Squares" using three different versions of Abel's name. Are there jet planes in a jet space? What was Voevodsky's motive; probably not fame. Do short Abelian squares make better company?

Some things I do understand which also have great names are the "classic" curves. I think the Quadratrix of Hippias sounds great and I want one to ride around on. Compared to its name it is a sightly disappointing shape. You can use it for trisecting angles and squaring cubes, but since it can't be constructed with ruler and compass it doesn't count as a solution, mathematicians are very picky. Hippias was, according to Plato, vain and boastful. Which goes to show, if you are vain and boastful, invent something useful so people will still be talking about you two and a half thousand years later but write your own biography or people will still mention how vain and boastful you are.

The Conchoid of Nicomedes takes its name from another Greek, who seems to have left only this to the world with no record of how vain and boastful he was. The word conchoid is supposed to recall a conch shell, the idea being that if you stack a lot of these they look something like the patterns on a mussel. Archimedes and his spiral are very famous compared to the other two but despite writing a whole book about them they were first studied by Conon (the barbarian). Often priority is not enough to get your name attached to something, even if then Spiral of Conon sounds better.

The second row are all examples of a family of Conchoids called Conchoids of de Sluze, named after a famous Belgian. The word Cissoid means ivy-like which it isn't. Strophoid means belt with a twist which is a strange thing to have a single word for. Trisectrix is the feminine curve that divides an angle in three, much like the Quadratrix of Hippias but a little more complicated.

The Witch of Agnesi is called a Lorentzian by physicists. Why it's called a Witch is apparently a (hilarious) pun in Latin and probably stuck because Agnesi was a woman and 18$^{th}$ century mathematicians weren't too sensitive to the gender bias in their field. Tschirnhausen, as well as having a very difficult name to spell, invented porcelain in Europe, which is quite a bit more impressive than this funny curve but unfortunately he didn't market himself well enough and someone else's name became attached to the discovery. Bernoulli (all of them) is very famous, the word lemniscate isn't but should be, it means decorated with ribbons.

Quadratrix of Hippias

Conchoid of Nicomedes

Spiral of Archimedes

Cissoid of Dioceles

Strophoid of Newton

Trisectrix of Maclaurin

Witch of Agnesi

Cubic of Tschirnhausen

Lemniscate of Bernoulli

A note on the construction of these for those who are interested (if you are really interested you can right click and view the source and see that I should have defined some more classes and functions to make my life easier when drawing these):

• Quadratrix of Hippias: Move a line segment down and a radius around a quadrant at a uniform rate so they reach the bottom at the same time. The intersection of radius and line defines the curve.
• Conchoid of Nicomedes: Draw a line from a point O to a point on a line X and add an extra distance k with the same slope as X varies down the line.
• Spiral of Archimedes: Starting from zero rotate around the origin and draw a line whose length is proportional to the angle. The end of the line traces out the spiral.
• Cissoid of Dioceles: Given a point O and a line L draw a line parallel to L through O. Let point on P vary on L and the projection of P onto the line through O be Q. R is the intersection of the line through Q perpendicular to OP and the points R trace the curve.
• Strophoid of Newton: The position of the orthocentre of a triangle with base given by a radius and the third vertex varying around the circle.
• Trisectrix of Maclaurin: The intersection points of two lines; one spinning three times as fast as the other.
• Witch of Agnesi: Draw a line from the base of the circle to the upper line and construct a right triangle with hypotenuse given by the line segment outside the circle, the third vertex traces the curve.
• Cubic of Tschirnhausen: The negative pedal curve of the parabola (look it up).
• Lemniscate of Bernoulli: This one is neat. I am using something called Watt's linkage. We stick three rigid rods together: a small one which rotates uniformly around a circle $a$; a longer one fixed at the centre of a different circle $b$ and a final even longer one $c$. We attach one end of $c$ to the free point of $b$ and the midpoint of $c$ to the free end of $a$. As $a$ rotates the free end of $c$ traces the curve!

Advanced String Theory

This entry is quite long; skip to the end to look at a piece of string if you don't feel like doing any maths. The Lagrangian is the difference between kinetic and potential energy terms. \[ L = T - V \] All of classical mechanics can be done using the Lagrangian and the equations of motion it obeys, it is completely equivalent to Newton's equations. If you have never seen a Lagrangian before you will have to take my word that it satisfies the equation \[ \frac{d}{dt} \frac{\partial L}{\partial \dot{\theta}_i } = \frac{\partial L}{\partial \theta_i } \] where the dot denotes a time derivative and the $\theta_i$ are the co-ordinates relevant for the problem. We're going to use pendulums which is why I chose the symbol $\theta$. Deriving this takes a while and involves concepts like D'Alembert's principle principle and variational calculus.

Accepting Lagrange's equations we can expand out the $\frac{d}{dt}$ to get, \[ \sum_j \frac{\partial^2 L}{\partial \dot{\theta}_i \partial \theta_j }\dot{\theta}_j + \sum_j \frac{\partial^2 L}{\partial \dot{\theta}_i \partial \dot{\theta}_j }\ddot{\theta}_j = \frac{\partial L}{\partial \theta_i }. \] The second partial derivatives have two indices, $i$ and $j$ and act on the one index object $\theta_j$ in the same way a matrix multiplies a vector. We define the two matrices \[ M_{ij} = \frac{\partial^2 L}{\partial \dot{\theta}_i \partial \dot{\theta}_j } \] the "mass" (for zero potential this would just be a diagonal matrix with the particle masses on the diagonal) and \[ F_{ij} = \frac{\partial^2 L}{\partial \dot{\theta}_i \partial \theta_j } \] The equation is now, \[ F \dot{\theta} + M \ddot{ \theta } = \frac{\partial L}{\partial \theta} \] I have dropped the indices and you should think of $\theta$ as a vector.

As in the previous post for numerical analysis we prefer first order equations, even if there are more of them. We define $\omega = \dot{\theta} $ and re-write the equation as, \[ \frac{d}{dt} \left( \begin{array}{c} \theta \\ \omega \end{array} \right) = \left( \begin{array}{cc} 0 & I \\ 0 & -M^{-1}F \end{array} \right) \left( \begin{array}{c} \theta \\ \omega \end{array} \right) + \left( \begin{array}{c} 0 \\ M^{-1} \frac{\partial L}{\partial \theta} \end{array} \right) \] As bad as it looks this is now in a form amenable to numerical analysis. We just have to calculate, $M$, $F$ and $\frac{\partial L}{\partial \theta}$. Let's get to it!

We'll make things a bit simpler for ourselves by making all the pendulums of equal length and all the pendulum bobs of equal mass. The following diagram summarises the situation for two pendulums.

From the diagram we can see (after staring at it for a while) that \[ \begin{array}{cc} x_1 = l \sin \theta_1 & y_1 = L \cos \theta_1 \\ x_2 = l (\sin \theta_1+\sin \theta_2) & y_2 = l (\cos \theta_1 + \cos \theta_2) \\ \vdots & \vdots \end{array} \] We will use the abbreviations \[ \begin{array}{c} c_i = \cos \theta_i \\ s_i = \sin \theta_i \end{array} \] The potential of any of the masses is $-mgh$ where $h$ is the height. This gives for all $N$ pendula a potential, \[ \begin{array}{l} V = -mg\left[ y_1 + y_2 + \ldots + y_N \right] \\ V = -mgl\left[ c_1 + (c_1 + c_2) + \ldots + (c_1 + \ldots + c_N) \right] \\ V = -mgl\left[ N c_1 + (N-1) c_2 + \ldots + c_N \right] \end{array} \] The kinetic energy is \[ T = \frac{m}{2} \left[ \dot{x}_1^2 + \dot{y}_1^2 + \dot{x}_2^2 + \dot{y}_2^2 + \ldots + \dot{x}_N^2 + \dot{y}_N^2 \right]. \] From the expressions for $x_i$ and $y_i$ we obtain, \[ \begin{array}{cc} \dot{x}_1 = l \cos \theta_1 \dot{\theta}_1 & \dot{y}_1 = l \sin \theta_1 \dot{\theta}_1\\ \dot{x}_2 = l (\cos \theta_1 \dot{\theta}_1 + \cos \theta_2 \dot{\theta}_2 ) & \dot{y}_2 = l (\sin \theta_1 \dot{\theta}_1 + \sin \theta_2 \dot{\theta}_2 ) \\ \vdots & \vdots \end{array} \] and the kinetic energy, \[ \begin{array}{l} T = \frac{m}{2} \left[ \dot{x}_1^2 + \dot{y}_1^2 + \dot{x}_2^2 + \dot{y}_2^2 + \ldots + \dot{x}_N^2 + \dot{y}_N^2 \right].\\ T = \frac{ml^2}{2} \left[ (c_1 \omega_1)^2 + (c_1 \omega_1 + c_2 \omega_2 )^2 + \ldots + (c_1 \omega_1 + c_2 \omega_2 + \ldots c_N \omega_N ) + c \leftrightarrow s \right]. \end{array} \]

This is a bit complicated so we will use matrix notation to keep track of things. With \[ J_{ab} = \begin{array}{ll} (N+1 - b) & \quad a < b \\ (N+1 - a) & \quad a \geq b \end{array} \] a symmetric matrix (write it out) we have, \[ T = \frac{ml^2}{2} \sum_{a = 1}^{N} \sum_{b = 1}^{N} J_{ab} c_a c_b \omega_a \omega_b+ c \leftrightarrow s = \sum_{a = 1}^{N} \sum_{b = 1}^{N} J_{ab} c_{ab} \omega_a \omega_b \] where \[ \begin{array}{c} c_{ab} = \cos (\theta_a - \theta_b) = c_{ba} \\ s_{ab} = \sin (\theta_a - \theta_b) = -s_{ba} \end{array} \]

Let's start by calculating $M$ and $F$. \[ \frac{ \partial T}{\partial \omega_j } = ml^2 \sum_{a = 1}^{N} \omega_a J_{aj} c_{aj} \] Then \[ M_{ij} = \frac{ \partial^2 T}{\partial \omega_i \partial \omega_j } = ml^2 J_{ij} c_{ij} \] and using the symetry of mixed partial derivatives \[ F_{ij} = \frac{ \partial^2 T}{\partial \theta_j \partial \omega_i } = ml^2 \left[ J_{ij} \omega_j s_{ij} - \delta_{ij} \sum_{a = 1}^{N} \omega_a J_{ai} s_{ia} \right] \] Now we need $\frac{\partial L}{\partial \theta_i}$. We have \[ \frac{\partial V}{\partial \theta_i} = mgl J_{1i} s_i \] and \[ \frac{\partial T}{\partial \theta_i} = ml^2 \sum_{a=1}^N J_{ai} s_{ai} \omega_a \omega_i \] Which gives $\frac{\partial L}{\partial \theta_i}$. Phew!

There actually is one more term. Adding damping in the Lagrangian framework is a little bit tricky, a full explanation is here. This gives us an extra parameter $\beta$ which we use to control the amount of drag. Below is a simulation of N coupled pendulua for $N = 1 \ldots 12$. Presumably you could let $N \rightarrow \infty$ and get a continuous string, I'll leave that to you.

N:

β: