guest post by Glyn Adgie and Tim van Beek
Or: Are the Milankovich Cycles Causing the Ice Ages?
The Milankovich cycles are periodic changes in how the Earth orbits the Sun. One question is: can these changes can be responsible for the ice ages? On the first sight it seems impossible, because the changes are simply too small. But it turns out that we can find a dynamical system where a small periodic external force is actually strengthened by random ‘noise’ in the system. This phenomenon has been dubbed ‘stochastic resonance’ and has been proposed as an explanation for the ice ages.
In this blog post we would like to provide an introduction to it. We will look at a bistable toy example. And you’ll even get to play with an online simulation that some members of the Azimuth Project created, namely Glyn Adgie, Allan Erskine and Jim Stuttard!
Equations with noise
In This Weeks Finds, back in ‘week308’, we had a look at a model with ‘noise’. A lot of systems can be described by ordinary differential equations:
If 
where the 
A way to make the above mathematically precise is to write the equation in the form
as an Ito stochastic differential equation driven by a Wiener process.
The mathematics behind this is somewhat sophisticated, but we don’t need to delve into it in this blog post. For those who are interested in it, we have gathered some material here:
• Stochastic differential equation, Azimuth Library.
Physicists and engineers who model systems with ordinary and partial differential equations should think of stochastic differential equations as another tool in the toolbox. With this tool, it is possible to easily model systems that exhibit new behavior.
What is stochastic resonance?
Listening to your friend who sits across the table, the noise around you may become a serious distraction, although humans are famous for their ability to filter the signal—what your friend is saying – from the noise—all other sounds, people’s chatter, dishes clattering etc. But could you imagine a situation where more noise makes it easier to detect the signal?
Picture a marble that sits in a double well potential:
This graphic is taken from
• Roberto Benzi, Stochastic resonance: from climate to biology.
This review paper is also a nice introduction to the topic.
This simple model can be applied to various situations. Since we think about the ice ages, we can interpret the two metastable states at the two minima of the potential in this way: the minimum to the left at -1 corresponds to the Earth in an ice age. The minimum to the right at 1 corresponds to the Earth at warmer temperatures, as we know it today.
Let us add a small periodic force to the model. This periodic force corresponds to different solar input due to periodic changes in Earth’s orbit: the Milankovich cycles. There are actually several different cycles with different periods, but in our toy model we will consider only one force with one period.
If the force itself is too small to get the marble over the well that is between the two potential minima, we will not observe any periodic change in the position of the marble. But if we add noise to the motion, it may happen that the force together with the noise succeed in getting the marble from one minimum to the other! This is more likely to happen when the force is pushing the right way. So, it should happen with roughly the same periodicity as the force: we should see a ‘quasi-periodic’ motion of the marble.
But if we increase the noise more and more, the influence of the external force will become unrecognizable, and we will not recognize any pattern at all.
There should be a noise strength somewhere in between that makes the Milankovich cycles as visible as possible from the Earth’s temperature record. In other words, if we think of the periodic force as a ‘signal’, there should be some amount of noise that maximizes the ‘signal-to-noise ratio’. Can we find out what it is?
In order to find out, let us make our model precise. Let’s take the time-independent potential to be
This potential has two local minima at 
The constant 
Including noise leads to the equation of motion, which is usually called a Langevin equation by physicists and chemists:
For more about what a Langevin equation is, see our article on stochastic differential equations.
This model has two free parameters, 
It is possible to write programs that generate simulations aka numerical approximations to solutions of stochastic differential equations. For those interested in how this works, there is an addendum at the bottom of the blog post that explains this. We can use such a program to model the three cases of small, high and optimal noise level.
In the following graphs, green is the external force, while red is the marble position. Here is a simulation with a low level of noise:
As you can see, within the time of the simulation there is no transition from the metastable state at 1 to the one at -1. That is, in this situation Earth stays in the warm state.
Here is a simulation with a high noise level:
The marble jumps wildly between the states. By inspecting the graph with your eyes only, you don’t see any pattern in it, do you?
And finally, here is a simulation where the noise level seems to be about right to make the signal visible via a quasi-periodic motion:
Sometimes the marble makes a transition in phase with the force, sometimes it does not, sometimes it has a phase lag. It can even happen that the marble makes a transition while the force acts in the opposite direction!
Some members of the Azimuth crew have created an online model that calculates simulations of the model, for different values of the involved constants. You can see it here:
• Stochastic resonance example, Azimuth Code Project.
We especially thank Allan Erskine and Jim Stuttard.
You can change the values using the sliders under the graphic and see what happens. You can also choose different ‘random seeds’, which means that the random numbers used in the simulation will be different.
Above, we asked if one could calculate for a fixed
There is no general definition that makes sense for all kinds of systems. For our model system, let’s assume that the expectation value of 
Then one useful definition of the signal to noise ratio 
So, can one calculate the value of 
If you would like to know more, have a look at this page:
• Stochastic resonance, Azimuth Library.
Also, read on for more about how we numerically solved our stochastic differential equation, and the design decisions we had to make to create a simulation that runs on your web browser.
For übernerds only
Numerically solving stochastic differential equations
Maybe you are asking yourself what is behind the online model? How does one ‘numerically solve’ a stochastic differential equation and in what sense are the graphs that you see there ‘close’ to the exact solution? No problem, we will explain it here!
First, you need C code like this:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define znew ((z=36969*(z&65535)+(z>>16)) << 16)
#define wnew ((w = 18000*(w&65535)+(w>>16))&65535)
static unsigned long z = 362436069, w = 521288629;
float rnor() {
float x, y, v;
x = ((signed long) (znew + wnew)) * 1.167239E-9;
if (fabs(x) < 1.17741)
return (x);
y = (znew + wnew) * 2.328306E-10;
if (log(y) < 0.6931472 - 0.5 * x * x)
return x;
x = (x > 0) ? 0.8857913 * (2.506628 - x) : -0.8857913 * (2.506628 + x);
if (log(1.8857913 - y) < 0.5718733 - 0.5 * (x * x))
return (x);
do {
v = ((signed long) (znew + wnew)) * 4.656613E-10;
x = -log(fabs(v)) * 0.3989423;
y = -log((znew + wnew)* 2.328306E-10);
} while (y + y < x * x);
return (v > 0 ? 2.506628 + x : 2.506628 - x);
}
int main(void) {
int index = 0;
for(index = 0; index < 10; index++) {
printf("%f\n", rnor());
}
return EXIT_SUCCESS;
}
Pretty scary, huh? Do you see what it does? We don’t think anyone stands a chance to understand all the ‘magical constants’ in it. We pasted it here so you can go around and show it to friends who claim that they know the programming language C and understand every program in seconds, and see what they can tell you about it. We’ll explain it below.
Since this section is for übernerds, we will assume that you know stochastic processes in continuous time. Brownian motion on 
![C[0, T].](https://s0.wp.com/latex.php?latex=C%5B0%2C+T%5D.&bg=ffffff&fg=333333&s=0 "C[0, T].")
Another way to think about it is to start with numerical approximations of ordinary differential equations. When you are handed such an equation of the form
with an initial condition 
with 
It is possible to extend some schemes from ordinary differential equations to stochastic differential equations that involve certain random variables in every step. The best textbook on the subject that I know is this one:
• Peter E. Kloeden and Eckhard Platen, Numerical Solution of Stochastic Differential Equations, Springer, Berlin, 2010. (ZMATH review.)
The extension of the Euler scheme is very simple, deceptively so! It is much less straightforward to determine the extensions of higher Runge–Kutta schemes, for example. You have been warned! The Euler scheme for stochastic differential equations of the form
is simply
where all the 
• George Marsaglia and Wai Wan Tsang, The Monty Python method for generating random variables, ACM Transactions on Mathematical Software 24 (1998), 341–350.
While an approximation to an ordinary differential equation approximates a function, an approximation to a stochastic differential equation approximates a probability distribution on 
A discrete approximation 
It is said to be of strong order 
As we see from the definition, a strong approximation needs to approximate the path of the original process ![[0, T].](https://s0.wp.com/latex.php?latex=%5B0%2C+T%5D.&bg=ffffff&fg=333333&s=0 "[0, T].")
What about weak convergence?
A discrete approximation 
where we assume that all appearing expectation values exist and are finite.
It is said to be of weak order
One speaks simply of ‘weak convergence’ when
So weak convergence means that the approximation needs to approximate the probability distribution of the original process 
Our simple online model implements the Euler scheme to calculate approximations to our model equation. So now you know it what sense it is indeed an approximation to the solution process!
Implementation details of the interactive online model
There appear to be three main approaches to implementing an interactive online model. It is assumed that all the main browsers understand Javascript.
1) Implement the model using a server-side programme that produces graph data in response to parameter settings, and use Javascript in the browser to plot the graphs and provide interaction inputs.
2) Pre-compute the graph data for various sets of parameter values, send these graphs to the browser, and use Javascript to select a graph and plot it. The computation of the graph data can be done on the designers machine, using their preferred language.
3) Implement the model entirely in Javascript.
Option 1) is only practical if one can run interactive programs on the server. While there is technology to do this, it is not something that most ISPs or hosting companies provide. You might be provided with Perl and PHP, but neither of these is a language of choice for numeric coding.
Option 2) has been tried. The problem with this approach is the large amount of pre-computed data required. For example, if we have 4 independent parameters, each of which takes on 10 values, then we need 10000 graphs. Also, it is not practical to provide much finer control of each parameter, as this increases the data size even more.
Option 3) looked doable. Javascript has some nice language features, such as closures, which are a good fit for implementing numerical methods for solving ODEs. The main question then is whether a typical Javascript engine can cope with the compute load.
One problem that arose with a pure Javascript implementation is the production of normal random deviates required for the numerical solution of an SDE. Javascript has a random() function, that produces uniform random deviates, and there are algorithms that could be implemented in Javascript to produce normal deviates from uniform deviates. The problems here are that there is no guarantee that all Javascript engines use the same PRNG, or that the PRNG is actually any good, and we cannot set a specific seed in Javascript. Even if the PRNG is adequate, we have the problem that different users will see different outputs for the same parameter settings. Also, reloading the page will reseed the PRNG from some unknown source, so a user will not be able to replicate a particular output.
The chosen solution was to pre-compute a sequence of random normal deviates, using a PRNG of known characteristics. This sequence is loaded from the server, and will naturally be the same for all users and sessions. This is the C++ code used to create the fixed Javascript array:
#include <boost/random/mersenne_twister.hpp>
#include <boost/random/normal_distribution.hpp>
#include <iostream>
int main(int argc, char * argv[])
{
boost::mt19937 engine(1234); // PRNG with seed = 1234
boost::normal_distribution<double> rng(0.0, 1.0); // Mean 0.0, standard deviation 1.0
std::cout << "normals = [\n";
int nSamples = 10001;
while(nSamples--)
{
std::cout << rng(engine) << ",\n";
}
std::cout << "];";
}
Deviates with a standard deviation other than 1.0 can be produced in Javascript by multiplying the pre-computed deviates by a suitable scaling factor.
The graphs only use 1000 samples from the Javascript array. To see the effects of different sample paths, the Javascript code selects one of ten slices from the pre-computed array.
Posted by John Baez 
measured in Celsius. We’ll assume the Earth radiates power square meter equal to
watts/meter2 and 
watts/meter2 per degree Celsius. This is a linear approximation taken from satellite data on our Earth. In reality, the power emitted grows faster than linearly with temperature.
is the average 
is the coalbedo: the fraction of solar power that gets absorbed. The coalbedo depends on the temperature; we’ll have to say how.
so that
but it’s easier to solve for 
The coalbedo, remember, is the fraction of sunlight that gets absorbed when it hits the Earth. It’s 1 minus the 
and 
Let’s try this:
this may seem scary. But don’t be intimidated!
is just a function that goes smoothly from 0 at low temperatures to 1 at high temperatures. This ensures that the coalbedo is near its icy value 
at high temperatures. But the fun part here is 
a parameter that says how rapidly the coalbedo rises as the Earth gets warmer. Depending on this, we’ll get different effects!
since that’s where 
has the biggest slope. We’re just lucky that in Celsius 
is the melting point of ice, so this makes a bit of sense.
In this graph the temperature ranges from -50 °C and 50 °C; the graph depends on what value of 
since then we get a graph like this:
on our slider. Then the coalbedo of the Earth becomes independent of temperature: it’s 0.525, halfway between its icy and ice-free values. Then, when the insolation takes its actual value of 342.5 watts per square meter, the model says the Earth’s temperature is very chilly: about -20 °C!
If we plug that in our formula:
were determined by satellite measurements of how much radiation actually escapes the Earth’s atmosphere and shoots out into space. As a further check to our model, we can look at an even simpler zero-dimensional energy balance model: a completely black Earth with no greenhouse effect. Another member of the Azimuth Project has written about this:
. When we start at some time 
on a given manifold 
, the flow of every fluid package is described by a path on 
there is a 
defined by the requirement that it maps the initial position 
at time 
. In order to apply geometric concepts in this situation, we will have to turn 
a neighborhood 
and an isomorphism (a “chart”) of 
, the “model space”. These isomorphisms can be used to transport concepts from 
should be taken in place of 
consisting of smooth maps from 
and 
open and 
a continuous map. The derivative of 
at the point 
in the direction 
is the map
![P: C^{\infty}[a, b] \to C^{\infty}[a, b]](https://s0.wp.com/latex.php?latex=P%3A+C%5E%7B%5Cinfty%7D%5Ba%2C+b%5D+%5Cto+C%5E%7B%5Cinfty%7D%5Ba%2C+b%5D+&bg=ffffff&fg=333333&s=0 "P: C^{\infty}[a, b] \to C^{\infty}[a, b]")
of a Fréchet manifold 
consists of all pairs 
where 
is a smooth curve
forms a Fréchet manifold, the tangent bundle 
, just as in finite dimensions.
is a smooth map
on 
:
:
of 
– to any other tangent space 
. If you can define an inner product on the Lie algebra, you can use this trick to transport the inner product to all the other tangent spaces by left or right multiplication, which will get you a left or right invariant metric.
of a tangent space 
there are unique vectors 
that are mapped to 
, that is
to have the value of that of 
:
, too, of course.
:
![\nabla_X Y - \nabla_Y X = [X, Y] \textrm{\; (torsion freeness)}](https://s0.wp.com/latex.php?latex=%5Cnabla_X+Y+-+%5Cnabla_Y+X+%3D+%5BX%2C+Y%5D+%5Ctextrm%7B%5C%3B+%28torsion+freeness%29%7D+&bg=ffffff&fg=333333&s=0 "\nabla_X Y - \nabla_Y X = [X, Y] \textrm{\; (torsion freeness)}")
![2 \langle \nabla_X Y, Z \rangle = X \langle Y, Z \rangle + Y \langle Z, X \rangle - Z \langle X, Y \rangle + \langle [X, Y], Z \rangle - \langle [Y, Z], X \rangle + \langle [Z, X], Y \rangle](https://s0.wp.com/latex.php?latex=2+%5Clangle+%5Cnabla_X+Y%2C+Z+%5Crangle+%3D+X+%5Clangle+Y%2C+Z+%5Crangle+%2B+Y+%5Clangle+Z%2C+X+%5Crangle+-+Z+%5Clangle+X%2C+Y+%5Crangle+%2B+%5Clangle+%5BX%2C+Y%5D%2C+Z+%5Crangle+-+%5Clangle+%5BY%2C+Z%5D%2C+X+%5Crangle+%2B+%5Clangle+%5BZ%2C+X%5D%2C+Y+%5Crangle+&bg=ffffff&fg=333333&s=0 "2 \langle \nabla_X Y, Z \rangle = X \langle Y, Z \rangle + Y \langle Z, X \rangle - Z \langle X, Y \rangle + \langle [X, Y], Z \rangle - \langle [Y, Z], X \rangle + \langle [Z, X], Y \rangle")
![2 \langle \nabla_X Y, Z \rangle = \langle [X, Y], Z \rangle - \langle [Y, Z], X \rangle + \langle [Z, X], Y \rangle](https://s0.wp.com/latex.php?latex=2+%5Clangle+%5Cnabla_X+Y%2C+Z+%5Crangle+%3D+%5Clangle+%5BX%2C+Y%5D%2C+Z+%5Crangle+-+%5Clangle+%5BY%2C+Z%5D%2C+X+%5Crangle+%2B+%5Clangle+%5BZ%2C+X%5D%2C+Y+%5Crangle+&bg=ffffff&fg=333333&s=0 "2 \langle \nabla_X Y, Z \rangle = \langle [X, Y], Z \rangle - \langle [Y, Z], X \rangle + \langle [Z, X], Y \rangle")
![\mathrm{ad}_X Z = [X, Z]](https://s0.wp.com/latex.php?latex=%5Cmathrm%7Bad%7D_X+Z+%3D+%5BX%2C+Z%5D+&bg=ffffff&fg=333333&s=0 "\mathrm{ad}_X Z = [X, Z]")
![\langle \mathrm{ad}^*_X Y, Z \rangle := \langle Y, \mathrm{ad}_X Z \rangle = \langle Y, [X, Z] \rangle](https://s0.wp.com/latex.php?latex=%5Clangle+%5Cmathrm%7Bad%7D%5E%2A_X+Y%2C+Z+%5Crangle+%3A%3D+%5Clangle+Y%2C+%5Cmathrm%7Bad%7D_X+Z+%5Crangle+%3D+%5Clangle+Y%2C+%5BX%2C+Z%5D+%5Crangle+&bg=ffffff&fg=333333&s=0 "\langle \mathrm{ad}^*_X Y, Z \rangle := \langle Y, \mathrm{ad}_X Z \rangle = \langle Y, [X, Z] \rangle")
![\mathrm{ad}_X Z = [X,Z]](https://s0.wp.com/latex.php?latex=%5Cmathrm%7Bad%7D_X+Z+%3D+%5BX%2CZ%5D&bg=ffffff&fg=333333&s=0 "\mathrm{ad}_X Z = [X,Z]")
. Another is to use ‘adjoint’ for the linear map 
coming from a linear map between inner product spaces 
given by 
Please don’t blame me for this.
and get
as an additional summand. So, in this case we get as the geodesic equation (again, for an invariant metric):
and its diffeomorphism group 
The Lie algebra 
of 
can be identified with the space of all vector fields on 
If we sloppily identify 
with coordinate 
and 
the commutator ![[X, Y] = (u v_x - u_x v) \partial_x](https://s0.wp.com/latex.php?latex=%5BX%2C+Y%5D+%3D+%28u+v_x+-+u_x+v%29+%5Cpartial_x+&bg=ffffff&fg=333333&s=0 "[X, Y] = (u v_x - u_x v) \partial_x")
is short for the derivative:
, this results in
for the involved functions.
here. In 2004, 
